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Abstract Now that detection of gravitational wave signals from the coales- 
cence of extra-galactic compact binary star mergers has become nearly routine, 
it is intriguing to consider other potential gravitational wave signatures. Here 
we examine the prospects for discovery of continuous gravitational waves from 
fast-spinning neutron stars in our own galaxy and from more exotic sources. 
Potential continuous-wave sources are reviewed, search methodologies and re- 
sults presented and prospects for imminent discovery discussed. 
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1 Introduction 


The LIGO and Virgo gravitational wave detectors have made historic discov- 
eries over the last seven years. The first direct detection in September 2015 
of gravitational waves marked a milestone in fundamental science (Abbott 
et al., 2016b), confirming a longstanding prediction of Einstein’s General The- 
ory of Relativity (Einstein, 1916, 1918). That the detection came from the 
first observation of a binary black hole merger provided a bonus not only in 
verifying detailed predictions of General Relativity, but in establishing unam- 
biguously that stellar-mass black holes exist in the Universe. More than 80 
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binary black hole (BBH) systems have been observed since GW150914 (Ab- 
bott et al., 2016d, 2017h,i,j, 2019b, 2021f,¢). Merging binary neutron star 
(BNS) systems (Abbott et al., 2017k, 2020a) have also been observed, includ- 
ing GW170817 (Abbott et al., 2017k), which was accompanied by a multitude 
of electromagnetic observations (Abbott et al., 20171). Those observations con- 
firmed the association of at least some short gamma ray bursts with binary 
neutron star mergers (Abbott ct al., 2017g) and the onset of kilonovae in BNS 
mergers that contribute substantially to the heavy element production in the 
Universe (Abbott et al., 20171). More recently came detections of merging neu- 
tron star — black hole (NSBH) systems (Abbott et al., 2021i). These discoveries 
of transient gravitational wave signals have ignited the field of gravitational 
wave astronomy. 


This review concerns a quite different and as-yet-undiscovered gravitational 
wave signal type, one defined by stability and near-monochromaticity over long 
time scales, namely continuous waves. CW signals with strengths detectable by 
current and imminent ground-based gravitational wave interferometers could 
originate from relatively nearby galactic sources, such as fast-spinning neutron 
stars exhibiting non-axisymmetry (Thorne, 1989), or more exotically, from 
strong extra-galactic sources, such as super-radiant Bose-Einstein clouds sur- 
rounding black holes (Arvanitaki et al., 2010). 


We already know from prior LIGO and Virgo searches that the strengths 
of CW signals must be exceedingly weak [~(10~7*) or less], which is consistent 
with theoretical expectation, from which we expect plausible CW strain am- 
plitudes to be orders of magnitudes lower than the amplitudes of the transient 
signals detected to date [~(10~21)]. This disparity in signal strength holds de- 
spite the much nearer distance of galactic neutron stars (~kpc) compared to 
the compact binary mergers (~40 Mpc to multi-Gpc) seen to date. In fact, it 
is only their long-lived nature that gives us any hope of detecting CW signals 
through integration over long data spans, so as to achieve a statistically viable 
signal-to-noise (SNR) ratio. As discussed below, however, that SNR increases, 
at best, as the square root of observation time, but for most CW searches, in- 
creases with an even lower power of observation time, while computational cost 
increases with much higher powers. These different scalings of signal sensitiv- 
ity and cost have led to a variety of approaches in targeting signals, depending 
on the size of signal parameter space searched. 


The search for continuous gravitational radiation has been under way since 
the 1970’s, using data from interferometers (Levine and Stebbins, 1972) and 
bars (Hirakawa et al., 1978; Suzuki, 1995), including from early prototypes (Li- 
vas, 1989) for the large gravitational wave detectors to come later. This review 
focuses primarily on the most recent searches from the Advanced LIGO and 
Virgo detectors, although summaries of search algorithm developments in the 
initial LIGO and Virgo era (and before) provide some historical context. For 
reference, the Advanced LIGO and Virgo runs to date comprise (with selected 
highlighted detections): 
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— The O1 observing run (LIGO only): September 12, 2015 — January 12, 2016 
— First detection of gravitational waves from a BBH merger: GW150914 (Ab- 
bott et al., 2016b). 

— The 02 observing run (LIGO joined by Virgo in last month): November 
30, 2016 — August 25, 2017 — First detection of gravitational waves from a 
BNS merger: GW170817 (Abbott ct al., 201 7k). 

— The O83 observing run (LIGO and Virgo): April 1, 2019 — March 27, 
2020 — First detection of gravitational waves from the formation of an 
intermediate-mass black hole: GW190521 (Abbott et al., 2020d) and the 
first detection of NSBH mergers. The run was divided into a 6-month “O3a” 
epoch (April 1, 2019 — October 1, 2019) and “O3b” (November 1 — March 
27, 2020) by a 1-month commissioning break. Many initial publications 
focused on results from the O38a data. 


In the following, section 2 reviews both conventional and exotic potential 
sources of CW gravitational radiation. Section 3 describes a wide variety of 
search methodologies being used to address the challenges of detection. Sec- 
tion 4 presents results (so far all only upper limits) from searches based on these 
algorithms, with an emphasis on the most recent results from the Advanced 
LIGO and Virgo detectors. Finally, section 5 discusses the outlook for discov- 
ery in the coming years, including the prospects for electromagnetic observa- 
tions of the continuous gravitational wave sources. This review focuses on CW 
radiation potentially detectable with current-generation and next-generation 
ground-based gravitational wave interferometers, which are sensitive to grav- 
itational frequencies in the human-audible band for sound. Past and future 
searches for lower-frequency CW radiation from supermassive black hole bi- 
naries at ~nHz frequencies using pulsar timing arrays (Manchester, 2010) or 
from stellar-mass galactic binaries at ~mHz frequencies using the space-based 
LISA (Bender et al., 1996) are not discussed here. 

Textbooks addressing gravitational waves, their detection and their analy- 
sis include references (Misner et al., 1972; Schutz, 1985; Maggiore, 2008, 2018; 
Saulson, 2017; Creighton and Anderson, 2011; Jaranowski and Krolak, 2009; 
Andersson, 2019). Review articles and volumes on gravitational wave science 
include references (Thorne, 1989; Blair et al., 1991; Sathyaprakash and Schutz, 
2009; Pitkin et al., 2011; Freise and Strain, 2010; Blair et al., 2012; Riles, 2013; 
Romano and Cornish, 2017). This review is a substantial expansion upon a 
briefer previous article (Riles, 2017). Other reviews of CW search methodol- 
ogy include (Prix, 2009; Palomba, 2012; Lasky, 2015; Sieniawska and Bejger, 
2019; Tenorio et al., 2021b; Piccinni, 2022). 


2 Potential sources of CW radiation 


In the frequency band of present ground-based detectors, the canonical sources 
of continuous gravitational waves are galactic, non-axisymmetric neutron stars 
spinning fast enough to produce gravitational waves in the LIGO and Virgo 
detectable band (at 1x, ~4/3x or 2x rotation frequency, depending on the 
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generation mechanism). These nearby neutron stars offer a “conventional” 
source of CW radiation — as astrophysically extreme as such objects are. 

A truly exotic postulated source is a “cloud” of bosons, such as QCD axions, 
surrounding a fast-spinning black hole, bosons that can condense in gargan- 
tuan numbers to a small number of discrete energy levels, enabling coherent 
gravitation wave emission from boson annihilation or from level transitions. 
Attention here focuses mainly on the conventional neutron stars, but the exotic 
boson cloud scenario is also discussed. 


2.1 Fast-spinning neutron stars 


The following subsections give an overview of neutron star formation, struc- 
ture, observables and populations, present the phenomenology of neutron-star 
spin-down, discuss potential sources of non-axisymmetry in neutron stars, and 
consider a number of particular GW search targets of interest. Although neu- 
tron stars were first postulated by Baade & Zwicky (Baade and Zwicky, 1934) 
and their basic properties worked out by Oppenheimer and Volkoff (Oppen- 
heimer and Volkoff, 1939), the first definitive establishment of their existence 
came with the discovery of the first radio pulsar (Hewish et al., 1968) PSR 
B1919+21 in 1967 with prior theoretical support for neutron star radiation 
contributing to supernova remnant shell energetics (Pacini, 1967) and rapid 
theoretical follow-up to explain the pulsation mechanism (Gold, 1968; Goldre- 
ich and Julian, 1969; Ruderman and Sutherland, 1975). 


2.1.1 Neutron star formation, structure, observables and populations 


As background, this section surveys at a basic level the fundamentals of neu- 
tron star formation, structure, observables and populations. Much more de- 
tailed information can be found in the following review articles or volumes on 
neutron stars (Lattimer and Prakash, 2001; Chamel and Haensel, 2008; Becker, 
2009; Ozel and Freire, 2016), pulsars (Lorimer and Kramer, 2005; Lyne and 
Graham-Smith, 2006; Lorimer, 2008), and rotating relativistic stars (Pascha- 
lidis and Stergioulas, 2017). 

Neutron stars are the final states of stars too massive to form white dwarfs 
upon collapse after fuel consumption and too light to form black holes, having 
progenitor masses in the approximate range 6-15 Mo (Lyne and Graham- 
Smith, 2006; Cerda-Duran and Elias-Rosa, 2018; Stockinger et al., 2020). 
These remarkably dense objects, supported by neutron degeneracy pressure, 
boast near-nuclear densities in their crusts and well-beyond-nuclear densities in 
their cores. The range of densities and associated total stellar masses and radii 
depend on an equation of state that is not experimentally accessible in terres- 
trial laboratories because of the combination of high density and (relatively) 
low temperature. A variety of equations of state have been proposed (Lat- 
timer and Prakash, 2001), with a small subset disfavored by the measurement 
of neutron stars greater than two solar masses (Buballa et al., 2014), by radii 
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of approximately ten kilometers (Miller et al., 2019b; Miller et al., 2021; Riley 
et al., 2019, 2021) and by the absence of severe tidal deformation effects in the 
gravitational waveforms measured for the BNS merger GW170817 (Abbott 
et al., 2017k, 2018c; Lim and Holt, 2019; Essick et al., 2020). The detection 
of a ~2.6-Mo object in the GW190814 merger (Abbott et al., 2020c) poses a 
challenge to the nuclear equation of state if the object is indeed a neutron star 
instead of a light black hole. 


In broad summary, a neutron star is thought to have a crust with outer 
radius between 10 and 15 km and a thickness of ~(1 km)(Shapiro and Teukol- 
sky, 1983), composed near the top of a tight lattice of neutron-rich heavy 
nuclei, permeated by neutron superfluid. Deeper in the star, as pressure and 
density increase, the nuclei may become distorted and elongated, forming a 
“nuclear pasta” of ordered nuclei and gaps (Ravenhall et al., 1983; Caplan and 
Horowitz, 2017). Still deeper, the pasta gives way to a hyperdense neutron 
fluid and perhaps undergoes phase transitions involving hyperons, perhaps to 
a quark-gluon plasma, or even perhaps to a solid strange-quark core (Shapiro 
and Teukolsky, 1983; Lattimer and Prakash, 2001). 


Uncertainties in equation of state lead directly to uncertainties in the ex- 
pected maximum mass and radius of a neutron star (Lattimer and Prakash, 
2001), but theoretical prejudice is consistent with the absence of observation in 
binary systems of neutron star masses much higher than two solar masses (Ozel 
and Freire, 2016; Demorest et al., 2010; Arzoumanian et al., 2018; Antoniadis 
et al., 2013; Cromartie et al., 2020; Fonseca et al., 2021). Neutron star radii 
are especially challenging to measure directly, with older measurements coming 
from X-ray measurements, where inferences are drawn from brightness of the 
radiation, its temperature and distance to the source, assuming black-body ra- 
diation, with corrections for the strong space-time curvature affecting the visi- 
ble surface area (Ozel and Freire, 2016; Degenaar and Suleimanov, 2018). New 
measurements from the NICER X-ray satellite are improving upon the preci- 
sion with which mass and radius can be determined simultaneously from indi- 
vidual stars, constraining more tightly the allowed equations of state (Miller 
et al., 2019b; Miller et al., 2021; Bogdanov et al., 2019a,b, 2021; Raaijmakers 
et al., 2019, 2021; Riley et al., 2019, 2021). Measurements of the gravitational 
waveform from the binary neutron star merger GW170817 have also provided 
new constraints and disfavor very stiff equations of state that lead to large neu- 
tron star radii (Abbott et al., 2018c). Detection of additional binary neutron 
star mergers in the coming years should improve these constraints. Broadly, 
one expects average neutron star densities of ~7 x 10!4 g cm~3, well above 
the density of nuclear matter (~3 x 1014) (Lorimer and Kramer, 2005), with 
densities at the core likely above 10'° g cm~3 (Shapiro and Teukolsky, 1983). 
See (Yunes et al., 2022) for a recent review of what has been learned about the 
neutron star equation of state from gravitational wave and X-ray observations. 
A recent Bayesian combined analysis (Huth et al., 2022) of predictions from 
chiral effective field theory of QCD, measured BNS gravitational waveforms, 
NICER X-ray observations and measurements from heavy ion (gold) collisions 
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indicate a somewhat stiffer equation of state than previously favored and hence 
larger allowed radii of neutron stars. 


Given the immense pressure on the nuclear matter, one expects a neu- 
tron star to assume a highly spherical shape in the limit of no rotation and, 
with rotation, to become an axisymmetric oblate spheroid. True axisymme- 
try would preclude emission of quadrupolar gravitational waves from rotation 
alone. Hence CW searchers count upon a small but detectable mass (or mass 
current) non-axisymmetry, discussed in detail in section 2.1.3. 


During the collapse of their slow-spinning stellar progenitors, neutron stars 
can acquire an impressive rotational speed as angular momentum conservation 
spins up the infalling matter. Even the slowest-rotating known pulsar spins on 
its axis every 24 seconds (Tan et al., 2018; Manchester and Hobbs, 2005), 
implying a rotational kinetic energy of ~(10°" J), and other young pulsars 
with spin frequencies of 10’s of Hz have rotational energies of ~(104? J). Recy- 
cled millisecond pulsars acquire even higher spins via accretion from a binary 
companion star, leading to measured spin frequencies above 700 Hz (Hessels 
et al., 2006; Bassa et al., 2017) and a rotational energy of ~(10*° J), or several 
percent of the magnitude of the gravitational bound energy of the star. This 
immense reservoir of rotational energy might appear to bode well for sup- 
porting detectable gravitational wave emission, but vast energy is required to 
create appreciable distortions in highly rigid space-time. From the perspective 
of gravitational-wave energy density (Misner et al., 1972), one can define an 


effective, frequency-dependent Young’s modulus Yom ~ Sow (~10%! Pa for 
few ® 100 Hz, or 20 orders of magnitude higher than steel). As a result, one 
must tap a significant fraction of the reservoir’s energy loss rate in order to 
produce detectable radiation, as quantified below, 


Most of the ~3300 known neutron stars in the galaxy are pulsars, detected 
via pulsed electromagnetic emission, primarily in the radio band, but also in X- 
rays and y-rays (with a small number detected optically) (Lyne and Graham- 
Smith, 2006; Manchester and Hobbs, 2005). Pulses are typically observed at 
the rotation frequency of the star, as a beam of radiation created by curvature 
radiation (Buschauer and Benford, 1976) from particles flung out from the 
magnetic poles (misaligned with the spin axis), sweeps across the Earth once 
per rotation (see (Melrose et al., 2021), however, for a critique of this model). 
A subset of neutron stars presumed to have magnetic poles tilted nearly 90 
degrees from the spin axis display two distinct pulses. 


Other neutron stars are known from detection of X-rays from thermal 
emission (heat from formation and perhaps from magnetic field decay), par- 
ticularly at sites consistent with the birth locations and times of supernova 
remnants (Lyne and Graham-Smith, 2006). Still other neutron stars are in- 
ferred from accretion X-rays observed in binary systems, particularly low-mass 
X-ray binaries with accretion disks (Lyne and Graham-Smith, 2006), although 
some accreting binaries with compact stars contain black holes, such as the 
high-mass X-ray binary Cygnus X-1. Figure 1 shows nearly the entire pop- 
ulation of currently known pulsars (Manchester and Hobbs, 2005) with spin 
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period P shorter than 20 seconds! in the P-P plane, where P is the first time 
derivative of the period. Red triangles show isolated pulsars, and blue circles 
show binary pulsars. 

Neutron stars have strong magnetic field intensities as a natural result of 
their collapse. If the magnetic flux is approximately conserved, the reduction 
of the outer surface of the star to a radius of ~(10 km) ensures a static surface 
field far higher than achievable in a terrestrial laboratory (Pacini, 1967), with 
inferred values (see below) ranging from 10° G to more than 10° G (Lyne and 
Graham-Smith, 2006). The strongest fields are seen in so-called “magnetars,” 
young neutron stars with extremly rapid spin-down, for which dynamo genera- 
tion is also likely relevant (Guilet and Miiller, 2015; Mésta et al., 2015). Both in 
young pulsars and in binary millisecond pulsars, there is reason to believe that 
stronger magnetic fields are “buried” in the star from accreting plasma (Payne 
and Melatos, 2004), although the burial mechanism is not confidently under- 
stood (Chevalier, 1989; Geppert et al., 1999; Lyne and Graham-Smith, 2006; 
Bernal et al., 2010). It has been suggested there is evidence in at least some pul- 
sars for slowly re-emerging (strengthening) magnetic field (Ho, 2011; Espinoza 
et al., 2011). Energy density deformation from a potentially non-axisymmetric 
buried field is another potential source of GW emission (Bonazzola and Gour- 
goulhon, 1996). See (Cruces et al., 2019) for a discussion of magnetic field 
decay preceding the accretion stage. 

In principle, there should be ~(108~°) neutron stars in our galaxy, (Narayan, 
1987; Treves et al., 2000). That only a small fraction have been detected is 
expected, for several reasons. Radio pulsations require high magnetic field and 
rotation frequency. Early studies (Goldreich and Julian, 1969; Sturrock, 1970; 
Ruderman and Sutherland, 1975; Lorimer and Kramer, 2005) implied the re- 
lation 

Bef = 1.7 X10" Ge (ig), (1) 


based on a model of radiation dominated by electron-positron pair creation in 
the stellar magnetosphere, a model broadly consistent with empirical observa- 
tion, although the resulting “death line” (see Figure 1) in the plane of period 
and period derivative is perhaps better understood to be a valley (Zhang et al., 
2000). 

The death line can be understood qualitatively from the following argu- 
ment. The rotating magnetic field of a neutron star creates a strong electric 
field that pulls charged particles from the star, forming a plasma with charge 
density po that satisfies (SI units): (Chen and Ruderman, 1993) 


po = —e0V - [(Q x r) x B(r)| (2) 
& —2e) Q-B(r), (3) 


where 2 is the angular velocity of the star, and B is the local magnetic field 
at location r with respect to the star’s center. In steady-state equilibrium, one 


1 The longest known pulsation period is 76 seconds from the recently discovered PSR 
J0901-4046 (Caleb et al., 2022), which also displays unusual pulse length and variability 
and which may represent a new pulsar class. 
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expects E- B = 0 since free charges can move along B-field lines. In so-called 
“gaps,” however, where the plasma density is low, a potential difference large 
enough to produce spontaneous electron-positron pair production can lead to 
synchroton RF radiation as the accelerated particles encounter curved mag- 
netic fields. This emission is thought to account for most radio pulsations (Lyne 
and Graham-Smith, 2006), where an “inner gap” refers to a region just outside 
the magnetic poles above the star’s surface, and an “outer gap” refers to a re- 
gion where a nominally dipolar magnetic field is approximately perpendicular 
to the rotation direction, separating regions of proton and electron flow from 
the star to the region beyond the “light cylinder,” defined by the cylindrical 
radius at which a co-rotating particle in the magnetosphere must travel at the 
speed of light. For the inner gap to have a voltage drop high enough to in- 
duce an amplifying cascade of pair production leading to coherent radio wave 
emission imposes a minimum value on the gap potential difference AV which, 
in general, can be approximated by (SI units): (Goldreich and Julian, 1969; 
Sturrock, 1970; Ruderman and Sutherland, 1975; Chen and Ruderman, 1993) 


(4) 


where F is the neutron star radius, leading (in a more detailed calculation) to 
equation 1 and via magnetic dipole emission assumptions (see section 2.1.2) 
to the death line shown in Figure 1 (but see (Smith et al., 2019) for evidence 
of selection effects and (Pétri, 2019) for a discussion of potentially important 
effects from higher order multipoles). Presumably, the vast majority of neu- 
tron stars created in the galaxy’s existence to date are now to the right of the 
line. Additional negative-sloped dashed lines in the figure indicate different 
nominal magnetic dipole field strengths and positive-sloped dashed lines indi- 
cate different nominal ages, based on observed present-day periods and period 
derivatives P/(2P) (see section 2.1.2). 

Two distinct major pulsar populations are apparent in Figure 1, defined 
by location in the diagram. The bulk of the population lies above and to the 
right of the line corresponding to B ~ 10! G. The bulk also lies above and to 
the left of the line corresponding to ages younger than ~(10°) years. Assuming 
a star’s magnetic field strength is stable, stars are expected to migrate down 
to the right along the B-field contours. Isolated pulsars seem to have typical 
pulsation lifetimes of ~(10" years), (Lyne and Graham-Smith, 2006) after 
which they become increasingly difficult to observe in radio. On this timescale, 
they also cool to where thermal X-ray emission is difficult to detect (Potekhin 
et al., 2015). There remains the possibility of X-ray emission from steady 
accretion of interstellar medium (ISM) (Ostriker et al., 1970; Blaes and Madau, 
1993), but it appears that the kick velocities from birth highly suppress such 
accretion (Hoyle and Lyttleton, 1939; Bondi and Hoyle, 1944) which depends 
on the inverse cube of the star’s velocity through the ISM, and steady-state X- 
ray emission from accretion onto even slow-moving neutron stars can be highly 
suppressed, consistent with non-observation to date of such accretion (Popov 
et al., 2015). 
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The remaining population, in the lower left of the figure, is characterized 
by shorter periods and smaller period derivatives. These are so-called “mil- 
lisecond pulsars” (MSPs), thought to arise from “recycling” of rotation speed 
due to accretion of matter from a binary companion. MSPs are stellar zom- 
bies, brought back from the dead with immense rotational energies imparted 
by infalling matter (Alpar et al., 1982; Radhakrishnan and Srinivasan, 1982). 
The rotation frequencies achievable through this spin-up are impressive — the 
fastest known rotator is PSR J1748—2446ad at 716 Hz (Hessels et al., 2006). 
One progenitor class for MSPs is the set of low mass X-ray binaries (LMXBs) 
in which the neutron star (~1.4 Mo) has a much lighter companion (~0.3 
Mo) (Lyne and Graham-Smith, 2006) that overfills its Roche lobe, spilling ma- 
terial onto an accretion disk surrounding the neutron star or possibly spilling 
material directly onto the star, near its magnetic polar caps. When the donor 
companion star eventually shrinks and decouples from the neutron star, the 
neutron star can retain a large fraction of its maximum angular momentum and 
rotational energy. Because the neutron star’s magnetic field decreases during 
accretion (through processes that are not well understood), the spin-down rate 
after decoupling can be very small. The minority of MSPs that are isolated are 
thought to have lost their one-time companions via consumption and ablation. 
A bridging class called “black widows” and “redbacks” refer to binary systems 
with actively ablating companions, such as B1957+20 (Fruchter et al., 1988; 
Strader et al., 2019; Roberts and van Leeuwen, 2013), where black widows 
denote the extreme subclass with companion masses below 0.1 Mg (Roberts 
and van Leeuwen, 2013). 


A nice confirmation of the link between LMXBs and recycled MSPs comes 
from “transitional millisecond pulsars” (tMSPs) in which accreting LMXB 
behavior alternates with detectable radio pulsations. The first tMSP found 
was PSR J1023+0038 (Bond et al., 2002; Thorstensen and Armstrong, 2005; 
Archibald et al., 2009), with two more systems since detected (Jaodand et al., 
2017). The nominal ages of MSPs extend beyond 10° years, that is, some 
have apparent ages greater than that of the galaxy (or even that of the Uni- 
verse). One possible explanation of this anomaly is reverse-torque spin-down 
during the Roche decoupling phase (Tauris, 2012), although a recent numeri- 
cal study suggests a more complex frequency evolution before and during the 
decoupling (Bhattacharyya, 2021). 


An obvious pattern in Figure 1, consistent with the recycling model, is 
the higher fraction of binary systems at lower periods. For example, binary 
systems account for 3/4 of the lowest 200 pulsar periods (below ~4 ms). 


Aside from the disappearance of stars from this diagram as they evolve 
toward the lower right and cease pulsations, there are also strong selection 
effects that suppress the visible population. We observe pulsars only if their 
radiation beams cross the Earth, only if that radiation is bright enough to 
be seen in the observing band, and only if the radiation is not sufficiently 
absorbed, scattered or frequency-dispersed to prevent detection with current 
radio telescopes. When the Square Kilometer Array project comes to fruition 
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Fig. 1 Measured rotational periods and period derivatives for known pulsars. Closed red 
triangles indicate isolated stars. Open blue circles indicate binary stars. The vertical dotted 
line denotes the approximate sensitivity band for Advanced LIGO at design sensitivity 
(faw > 10 Hz, assuming faw = 2frot). 


in the late 2020’s, it is estimated that the current known population of pulsars 
will grow tenfold (Kramer and Stappers, 2015). 


2.1.2 Neutron star spin-down phenomenology and mechanisms 


Nearly every known pulsar is observed to be spinning down, that is, to have 
a negative rotational frequency time derivative, implying loss of rotational 
kinetic energy. As discussed below in detail, there are many physical mech- 
anisms, electromagnetic and gravitational, that can lead to this energy loss. 
For CW signal detection we want a gravitational wave component, but there 
is good reason to believe that electromagnetic processes dominate for nearly 
every known pulsar. 

A convenient and commonly used phenomenological model for spin-down 
is a power law: 


fan (5) 
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where f is the star’s instantaneous frequency (rotational f;o¢ or gravitational: 
faw & frot), f is the first time derivative, and K is a negative constant for all 
but a handful of stars (thought to be experiencing large acceleration toward 
us because of nearness to a deep gravitational well, such as in the core of a 
globular cluster). The exponent n depends on the spin-down mechanism and is 
known as the braking index. The four most common theoretical braking indices 


discussed in the literature are the following: 


— n=1- “Pulsar wind” (extreme model) 

— n=3- Magnetic dipole radiation 

— n=5-— Gravitational mass quadrupole radiation (“mountain”) 

— n=7- Gravitational mass current quadrupole radiation (r-modes). 


In principle, other oscillation modes that can generate gravitational waves are 
also possible, but the n =5 and n=7 modes discussed below are thought to 
be the most promising. 

Assuming the same power law has applied since the birth of the star, the 
age 7 of the star can be related to its birth rotation frequency fp and current 
frequency f by (n 4 1): 


pe lee, f i)" e) 


and in the case that f < fo, 


A common baseline assumption in radio pulsar astronomy is that the brak- 
ing index is n = 3 from which the nominal magnetic dipole age of a star can 
be defined 


Tmag =r (8) 


again, under the assumption f < fo. 
From the more generic power-law spin-down model (Equation 5), the 2nd 
frequency derivative can be written: 


(ap far, (9) 


from which the current braking index can be determined if the spin frequency’s 
2nd time derivative can be measured reliably: 


fF 
Pe 
Before examining the empirical measurements of the braking indices, which 
are mostly inconsistent with n = 8, let’s briefly review spin-down mechanisms 
with well defined braking indices, when dominant. For GW radiation spin- 
down dominance, related “spin-down” limits on strain amplitude will also be 
presented. 


(10) 
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2.1.2.1 “Pulsar wind” (n = 1). 


Early on in pulsar astronomy (Michel, 1969; Michel and Tucker, 1969) 
it was recognized that the streaming of relativistic particles (electrons and 
positrons mainly, with some ions) away from the magnetosphere of a fast- 
spinning neutron star would lead to a spin-down torque that could, in principle, 
rival that from magnetic dipole radiation, in addition to distorting the shape 
of the magnetic field lines and affecting the dipole radiation (Gaensler and 
Slane, 2006). In this perhaps too-simple model, the spin-down is dominated 
by a braking torque from a return current (predominantly counter-flowing elec- 
trons and positrons) crossing magnetic field lines in the polar cap regions of 
the star (Contopoulos et al., 1999), leading to a braking index of one. A more 
recent study of magnetar spin-down (Harding et al., 1999) considered a model 
with sporadic high winds following bursts, with magnetic dipole emission dom- 
inating spin-down between bursts. In the steady state, however, considering 
the interaction of the magnetic field and the plasma of the magenetosphere, 
both magnetic dipole emission and pulsar wind contributions tend to yield 
a braking index of about three (Michel and Li, 1999; Spitkovsky, 2004), dis- 
cussed next. A phenomenological model (Melatos, 1997) that is a variant of 
the vacuum dipole mode, featuring an inner magnetosphere strongly coupled 
to the star, accounts successfully for the braking indices of the Crab and other 
young pulsars with n < 1. 


2.1.2.2 Magnetic dipole (n = 3). 


Equating rotational energy loss rate to magnetic dipole radiation losses, 
leads to the relation (Pacini, 1968): 


dE poM? w4 
Ean OTE ay) 


where w is the rotational angular speed, M, is the component of the star’s 
magnetic dipole moment perpendicular to the rotation axis (taken to be z 
axis): M, = Msin(a), with a the angle between the axis and north magnetic 
pole. In a pure dipole moment model, the magnetic pole field strength at 
the surface is Bo = oM /27R?. Equating this energy loss to that of the 
(Newtonian) rotational energy }J,,,w" leads to the prediction: 


du) oR? 9 3 
a “je. (12) 


Hence the magnetic dipole spin-down rate is proportional to the square of 
B, = Bosin(q) and to the cube of the rotation frequency, giving n = 3. 
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2.1.2.3 Gravitational mass quadrupole (“mountain”, n = 5). 


Let’s now consider the gravitational radiation one might expect from these 
stars. It is conventional to characterize a star’s mass quadrupole asymmetry 
by its equatorial ellipticity: 


a Ss Se, (13) 


An oblate spheroid naturally has a polar ellipticity, but in the absence of 
precession’, such a deformation does not lead to GW emission Henceforth 
“ellipticity” will refer to equatorial ellipticity, often attributed to a “mountain”. 
For a star at a distance d away and spinning about the approximate symmetry 
axis of rotation (z), (assumed optimal — pointing toward the Earth), then the 
expected intrinsic strain amplitude hg is 


2 2 
ho = 4m Gel fow (14) 


cr 
= in 10-°9 (a3) (32) (iit) (APE). os 


where Ip = 10°° kg-m? (104° g-cm?) is a nominal moment of inertia of a neutron 
star, and the gravitational radiation is emitted at frequency faw = 2 frot. The 
total power emission in gravitational waves from the star (integrated over all 
angles) is 


ae é te" (16) 


=r i019) (2) (Gs) (BR). an 


Equating this loss to the reduction of rotational kinetic energy 51 w?,, leads 
to the spin-down relation: 
; 3274 G ous 
iw == 5 Ps few (18) 
2 ( faw \” 
= (1.7 x 10-® Hz/s) (——) 1 
eet Ba) pe) toe * (19) 


in which the braking index of 5 is apparent. 
For an observed neutron star of measured f and f, one can define the “spin- 
down limit” on maximum allowed strain amplitude by equating the power loss 


2 Free precession of an oblate neutron star can lead to gravitational radiation at the 
rotation frequency (Zimmermann and Szedenits, 1979), but there is little empirical evidence 
for such precession in pulsars and good reason to expect that such precession would be 
rapidly damped by internal dissipation (Jones and Andersson, 2002). 
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in equation (16) to the time derivative of the (Newtonian) rotational kinetic 
energy: $1 w”, as above for magnetic dipole radiation. One finds: 


1 | 8G. fox 
Aspin— own — Loz 
spin—d d Ie faa 
_os, (1 kpe 1kHz —faw Lon 
= (2.6 x 107° 
pane 1( d ) (Gan) Gotan) (F2)| 


Hence for each observed pulsar with a measured frequency spin-down and 
distance d, one can determine whether or not energy conservation even per- 
mits detection of gravitational waves in an optimistic scenario. Unfortunately, 
nearly all known pulsars have strain spin-down limits below what can be de- 
tected by the LIGO and Virgo detectors at current sensitivities, as detailed 
below. 


Nie 
NO 
S 
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2.1.2.4 Gravitational mass current quadrupole (r-modes, n = 7). 


Different frequency scalings apply to mass quadrupole and mass current 
quadrupole emission. The most promising source of the mass current non- 
axisymmetry in neutron stars is thought to be “r-modes,” due to fluid motion 
of neutrons (or protons) in the crust or core of the star. Like jet streams in the 
Earth’s atmosphere that manifest Rossby waves, these currents are deflected 
by Coriolis forces, giving rise to spatial oscillations (Andersson, 1998; Bild- 
sten, 1998; Friedman and Morsink, 1998; Owen et al., 1998). These r-modes 
can be inherently unstable, arising from azimuthal interior currents that are 
retrograde in the star’s rotating frame, but which are prograde in an exter- 
nal reference frame. As a result, the quadrupolar gravitational wave emission 
due to these currents leads to an increase in the strength of the current. This 
positive-feedback loop leads to a potential intrinsic (Chandrasekhar-Friedman- 
Schutz (Chandrasekhar, 1970; Friedman and Schutz, 1978)) instability. The 
frequency of such emission is expected to be a bit more than approximately 
4/3 the rotation frequency (Andersson, 1998; Bildsten, 1998; Friedman and 
Morsink, 1998; Owen et al., 1998; Kojima, 1998; Caride et al., 2019). 

Following the notation of Owen (Owen, 2010; Caride et al., 2019), the mass 
current can be treated as due to a velocity field perturbation 6v,, integration 
over which leads to the following expression for the intrinsic strain amplitude 
seen at a distance d: 


[51207 G1. : 


_26 ( Lkpe few \?/ a R : 
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where a is the dimensionless r-mode amplitude, M is the stellar mass, RF its 
radius, and J is a dimensionless functional of the stellar equation of state, 
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which for a Newtonian polytrope with index 1 gives J % .0164 (Owen, 2010), 
assumed in the fiducial Eqn. 22. 
The energy loss in this model is (Thorne, 1980; Owen, 2010) 


dE 102479 G 


an ae fesxyoP M7 REF, (23) 
Equating this loss to the reduction of rotational kinetic energy ae as 
above, leads to the spin-down relation: 
, 4096 x’ G M?R°J? 
few = a” few (24) 


225 ch Ly 


~9.0 x 10-4 Hz/s ( —* "( c y fow_\" (25) 
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in which the braking index of 7 is apparent. 
As before, one can define a spin-down limit, but one based on pure r-mode 


radiation: 
1} 45-4 fay 
hs in—down — Liz ; 26 
F a r 8 ou foaw ( ) 


where the ratio of this spin-down limit to the one given in Equation 20 is 3/2, 
which arises simply from the different ratios of GW signal frequency to spin 
frequency for mass quadrupole vs. mass current quadrupole radiation®. 


I 


2.1.2.5 Measured braking indices. 


Figure 2 shows the distribution of 12 reliably measured braking indices from 
a recent snapshot of the ~3300 pulsars listed in the ATNF catalog (release 
V1.66 — January 10, 2022 (Manchester and Hobbs, 2005)). Nearly all have 
values below the nominal value of 3 for a magnetic dipole radiator, although 
several have large uncertainties. 

This distribution suggests that the model of a neutron star spinning down 
with constant magnetic field is, most often, inaccurate (Lyne and Graham- 
Smith, 2006). All measured values for this collection lie below 3.0, except X- 
ray pulsar PSR J1640—4631 with a measured index of 3.15+0.03 (Archibald 
et al., 2016). It is possible that for many stars the departure of the measured 
braking index from the nominal value is due to an admixture of magnetic 
dipole radiation and other steady-state processes (Melatos, 1997), although 
secular mechanisms may also play a role. See (Palomba, 2000, 2005) for dis- 
cussions of spin-down evolution in the presence of both gravitational wave and 
electromagnetic torques. Other suggested mechanisms for less-than-3 braking 
indices are decaying magnetic fields (Romani, 1990), re-emerging buried mag- 
netic fields (Ho, 2011), a changing inclination angle between the magnetic 


3 The 4/3 ratio assumed here for faw/ frot is a slow-rotation approximation in Newtonian 
gravity; the ratio changes by tens of percent for fast rotation in General Relativity (Idrisy 
et al., 2015; Caride et al., 2019) (see section 4.2). 
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Fig. 2 Measured braking indices inferred from frequency derivatives of young pulsars with 
rotation frequencies greater than 10 Hz. For frequently glitching pulsars, such as Vela, the 
braking index is computed as a long-term average (Espinoza et al., 2017). Horizontal bars 
indicate uncertainties and are smaller than the plot markers for several pulsars. Vertical 
lines at braking indices of 3 and 5 denote the nominal expectations for magnetic dipole and 
gravitational quadrupole emission, respectively. References: 1 (Lyne et al., 2015), 2 (Ferdman 
et al., 2015), 3 (Espinoza et al., 2017), 4 (Weltevrede et al., 2011), 5 (Clark et al., 2016), 
6 (Livingstone and Kaspi, 2011), 7 (Archibald et al., 2016), 8 (Espinoza et al., 2011), 9 (Roy 
et al., 2012), 10 (Livingstone et al., 2007). 


dipole axis the spin axis (Middleditch et al., 2006; Tauris and Konar, 2001; 
Ho, 2015; Lyne et al., 2015; Johnston and Karastergiou, 2017), and a changing 
superfluid moment of inertia (Ho and Andersson, 2012). An interesting obser- 
vation of the aftermath of two short GRBs noted indirectly inferred braking 
indices near or equal to three (Lasky et al., 2017a), suggesting the rapid spin- 
down of millisecond magnetars, possibly born from neutron star mergers. (No 
direct gravitational wave evidence of a such a post-merger remnant has been 
observed from GW170817 (Abbott et al., 20170, 2019d).) Similarly, a recent 
analysis of X-ray afterglows of gamma-ray bursts (Sarin et al., 2020) argues 
that at least some have millisecond magnetar remnants powering their emis- 
sion, with GRB 061121 yielding a braking index n = ABST Te consistent with 


18 Keith Riles 


JO857—4424 
J1841—0425 
J1824—1945 
J1320—5359 
J1806—2125 
J1715—3903 
J1028—5819 
J1731—4744 
J1718—3825 
J1420—6048 
JO835—4510 
J1531—5610 
J1524—5625 
J1112—6103 
J1357—6429 
J1801—2451 
J1709—4429 
J1826—1334 
J1637—4642 
J1648—4611 
J1617—5055 
J1803—2137 
J1048—5832 
J1830—1059 
J0940—5428 
J1646—4346 
J1301—6305 
J1809—1917 
J1016—5857 
JO908—4913 
J1410—6132 
J1730—3350 
J1412—6145 
J1726—3530 
J0954—5430 
J1015—5719 
J1643—4505 
J1737—3137 
J1614—5048 
J1702—4310 
JO901—4624 
J0659+1414 
J1702—4128 
J1509—5850 
J1815—1738 
J0855—4644 
J1632—4818 © ) Glitched 

J1524—5706 ~~ This work 
J1513—5908 — PJS+20 

J1734—3333 ° 


1 10 100 1000 
Braking index (n) 


Fig. 3 Measured braking indices inferred from frequency derivatives of pulsars compiled in 
(Lower et al., 2021) from (Parthasarathy et al., 2020) (PJS), and (Lower et al., 2021). An 
ensemble of glitching (open circles) and non-glitching pulsars are included. For most of the 
glitching stars, the braking indices are representative of their average inter-glitch braking, 
not their long-term evolution. Figure is taken from (Lower et al., 2021). 
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gravitational radiation dominance (albeit with large required ellipticity (Ho, 
2016; Kashiyama et al., 2016)). See (Strang et al., 2021), however, for an al- 
ternative study in which radiation driven from a millisecond magnetar can 
account for short GRB X-ray afterglows. 

It has been argued that the inter-glitch evolution of spin for the X-ray 
pulsar PSR J0537—6910 displays behavior consistent with a braking index 
of 7, (Andersson et al., 2018; Ho et al., 2020) consistent with r-mode emis- 
sion, while the long-term trends points to an underlying braking index of 
—1.25+0.01 (Ho et al., 2020). When intepreting the generally low values of well 
measured braking indices, one must bear in mind the potential for selection 
bias. Baysesian analysis of the spin evolution of 19 young pulsars (Parthasarathy 
et al., 2019, 2020), taking into account timing noise and extracting the long- 
term behavior from short-term, glitch-driven fluctuations, leads to braking 
indices much larger than 3. A similar follow-up analysis of an ensemble of 
glitching and non-glitching pulsars (Lower et al., 2021) confirmed that brak- 
ing indices exceeding 100 are observed (see Figure 3), albeit for stars in which 
a simple power-law spin-down is clearly inappropriate. 


2.1.2.6 The gravitar model and associated figures of merit. 


Gravitars refer to neutron stars with spin-down dominated by gravitational 
wave energy loss. Although there is good reason to believe that most known 
pulsars are not gravitars, nonetheless the model is useful in bounding expec- 
tation on what is possibly detectable. Figure 4 shows a subset of the pulsars 
from Figure 1, now graphed in the faw—faw plane, under the assumption that 
faw = 2 frot. Again, isolated and binary stars are denoted by closed circles 
and open triangles, respectively. A vertical dashed line bounds the approxi- 
mate detection bandwidth for Advanced LIGO at design sensitivity (~10 Hz 
and above). As in Figure 1, contours are shown for constant magnetic field, 
assuming spin-down dominated by magnetic dipole emission (n = 3). In addi- 
tion, contours of higher slope are shown for constant ellipticity An intriguing 
deficit of millisecond pulsars with extremely low period derivatives appears 
consistent (Woan et al., 2018) with a population of sources with a minimum 
ellipticity of about ~10~° with additional spin-down losses from magnetic 
dipole radiation (see near absence of sources in figure 4 to the right of the 
« = 10-® line). At the other extreme are lower-frequency, younger pulsars 
with high spin-downs, the highest of which is 7.6 x 107° Hz/s (Crab pulsar). 

Using Equation 20, these known pulsars can be mapped onto a plane of 
faw-—ho under the gravitar assumption, indicated in Figure 5. That is, the spin- 
down strain limit (for n = 5) is shown on the vertical axis. Also shown are 
corresponding contours of constant implied values of e/d, under the gravitar 
assumption, where d is the distance to the star. In addition, detector net- 
work sensitivities are shown for Advanced LIGO at design sensitivity (Abbott 
et al., 2020b) (two detectors for two observing years), henceforth designated 
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Fig. 4 Nominal expected GW frequencies and frequency derivatives for known pulsars. 
Closed triangles indicate isolated stars. Open circles indicate binary stars. Contours are 
shown for constant magnetic fields (ellipticities) for spin-down dominated by magnetic dipole 
(gravitational mass quadrupole) emissions. In this figure and in Figs. 5-8, the frequency 
derivatives have been corrected for the Shklovskii effect (Shklovskii, 1970) (apparent negative 
frequency derivative due to proper motion orthogonal to the line of sight). 


as the “O4/05 run” +, for two proposed configurations of the “3rd-generation” 
Einstein Telescope (Maggiore et al., 2020) (B and C, for three detectors for 
five observing years)” These sensitivities assumed a targeted search (discussed 
below) using known pulsar ephemerides. If a star is marked above a sensitiv- 
ity curve, then it is at least possible to detect it if its spin-down makes it a 


4 Although the O4 run scheduled to start near the start of 2023 will likely run for only 
~1 year (Abbott et al., 2020b) and may not quite reach the original Advanced LIGO design 
sensitivity, the succeeding O5 run in the “A+” configuration is expected to exceed Advanced 
LIGO sensitivity significantly and to last for more than a year, making the detector sensi- 
tivities assumed here conservative. On the other hand, the 04/05 observing time assumed 
here does not account for realistic deadtime losses. 

5 Another 3rd-generation proposal is for the “Cosmic Explorer” (Abbott et al., 2017c) 
which would have performance comparable to that of the Einstein Telescope (ET), being 
more sensitive at frequencies above ~10 Hz and less sensitive at lower frequencies. To avoid 
clutter in these figures, only the ET sensitivities are shown. 
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Fig. 5 Nominal expected GW frequencies and nominal strain spin-down limits for known 
pulsars. Closed triangles indicate isolated stars. Open circles indicate binary stars. The solid 
curves indicate the nominal (idealized) strain noise sensitivity for the O3 observing run 
(black), and expected sensitivities for 2-year Advanced LIGO data run at design sensitivity 
(magenta) and a 5-year Einstein Telescope data run for two different detector designs: 
ETB (blue) and ETC (green). Dashed diagonal lines correspond to particular quotients 
of ellipticity over distance. A subset of pulsars of particular interest are labeled on the 
figure. 


gravitar. Note, however, that Equation 20 has been applied with a nominal 
moment of inertia [,,, but the uncertainty in [,, if of order a factor of two, 
depending on equation of state and stellar mass (Worley et al., 2008). 
Another take on the pulsars with accessible spin-down limits is shown in 
Figure 6, where accessible ellipticity € values are shown for Advanced LIGO 
and Einstein Telescope (ETC) sensitivities. Each vertical bar represents a 
range of ellipticities detectable for that star (red = accessible to Advanced 
LIGO, green = accessible to Einstein Telescope), where the asterisk at the top 
of the each bar is the ellipticity corresponding to that star’s spin-down limit, 
given its few, fow and distance d values, while the depth to which the bar 
falls indicates the lowest detectable ellipticity. Straight dashed lines of nega- 
tive slope depict corresponding fow values under the mass quadrupole gravitar 
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Fig. 6 Nominal expected GW frequencies and maximum allowed ellipticities for known 
pulsars. Black or blue asterisks indicate ellipticities accessible with Advanced LIGO (2 de- 
tectors, 2 years) or ETC sensitivities (3 detectors, 5 years), respectively, using targeted 
searches, where red vertical lines terminated by red asterisks indicate ellipticity sensitivity 
range for Advanced LIGO, and green vertical lines and green asterisks indicate additional 
ellipticity sensitivity range for ETC. A selection of pulsars accessible with Advanced LIGO 
sensitivity are labeled in red. Diagonal dashed lines correspond to corresponding foaw values 
under the gravitar model. 


model. The actual faw may be significantly higher because of the spin-down 
mechanisms discussed earlier. A striking feature of this figure is that sensitiv- 
ities to very low ellipticities come almost entirely from the highest-frequency 
stars (as a reminder from Equation 14, ho « ef€y). For example, no known 
pulsar with an ellipticity below 107° and that is accessible to Advanced LIGO 
has a faw value lower than 70 Hz, and no ellipticity below 1078 is accessible 
to Advanced LIGO below 300 Hz. 

Another figure of merit is the distance to which searches can detect sources 
of a particular ellipticity. Figure 7 shows the estimated distances to known 
pulsars over the detection frequency band. Also shown are solid contours of 
Advanced LIGO sensitivity range for different ellipticity values and dashed 
contours for Einstein Telescope. Pulsars with spin-down limits accessible to 
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Fig. 7 Maximum allowed targeted-search ranges for gravitars vs GW frequencies for dif- 
ferent assumed ellipticities for Advanced LIGO sensitivity (solid magenta curves) and cor- 
responding ranges for ETC sensitivity (dashed blue curves). Known pulsar distances are 
shown vs the expected GW frequencies, where red dots indicate pulsars with accessible 
spin-down limits for Advanced LIGO sensitivity, and smaller green dots indicated pulsars 
with accessible spin-down limits for ETC sensitivity. Known pulsars in distinct horizontal 
bands (common distance) arise from stars in clusters or from distance capping in the galactic 
electron density model (Yao ct al., 2017) used in the ATNF catalog (Manchester and Hobbs, 
2005). 


Advanced LIGO are shown in red, and those accessible to Einstein Telescope 
are shown in green. Only a handful of pulsars within 500 pc are accessible to 
Advanced LIGO with ellipticities below 10~°. On the other hand, to reach the 
galactic center (~8.5 kpc) at a signal frequency of 1 kHz requires an ellipticity 
larger than ~3 x 1078, and at 100 Hz requires an ellipticity greater than 
~3 x 1078. 

As discussed in detail below, all-sky searches for unknown neutron stars 
necessarily have reduced sensitivity, such that the ranges shown for targeted 
searches using known pulsar timing do not apply. Figure 8 shows another range 
vs frequency plot, but for which (optimistic) Advanced LIGO and Einstein 
Telescope all-sky sensitivities are assumed. For reference, the all-sky strain 
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sensitivity is taken to be about 20 times worse than its targeted-search sensi- 
tivity for Advanced LIGO and the corresponding ratio about 40 times worse 
for Einstein Telescope. Consequently, the all-sky range contours corresponding 
to those in Figure 7 would be reduced by the same ratios. Alternatively, to 
obtain the same ranges in the all-sky search would require ellipticities higher 
by the same ratios. The all-sky ranges in Figure 8, in contrast, are shown 
as contours for different assumed faw values under the gravitar assumption. 
These contours are useful in assessing all-sky searches, since those searches 
are defined, in part, by their maximum spin-down range, which affects compu- 
tational cost. Once again, known pulsars for which this search technique can 
reach the spin-down limit are shown in red for Advanced LIGO and in green 
for Einstein Telescope. We see that for Advanced LIGO to reach the galactic 
center at a signal frequency of 1 kHz requires a minimum spin-down magnitude 
greater than 10~° Hz/s (minimum because another mechanism, such as mag- 
netic dipole emission, may contribute to a higher spin-down magnitude), and 
at 100 Hz requires a minimum spin-down magnitude just less than 10~!° Hz/s. 
The corresponding required ellipticities at those frequencies are ~8 x 10~° and 
~8 x 1077, respectively. 


A simple steady-state argument by Blandford (Thorne, 1989) led to an 
early estimate of the maximum detectable strain amplitude expected from a 
population of isolated gravitars of a few times 10-24, independent of typical 
ellipticity values, in the optimistic scenario that most neutron stars become 
gravitars. A later detailed numerical simulation (IKnispel and Allen, 2008) re- 
vealed, however, that the steady-state assumption does not generally hold for 
mass quadrupole radiation, leading to ellipticity-dependent expected maxi- 
mum amplitudes that can be 2-3 orders of magnitude lower in the LIGO band 
for ellipticities as low as 107° and a few times lower for ellipticity of about 
10~-°. Mass current quadrupole (r-mode) emission, however, would spin stars 
down faster, leading back to more optimistic maximum amplitudes (Owen, 
2010). A more detailed simulation including both electromagnetic and gravia- 
tional wave spin-down demonstrated the potential for setting joint constraints 
on natal neutron star magnetic fields and ellipticities (Wade et al., 2012). A 
recent population simulation study (Reed et al., 2021) estimated fractions of 
neutron stars probed by previous CW searches for different assumed elliptic- 
ities and concluded that the greatest potential gain from improving detector 
sensitity in accessing more neutron stars of plausible ellipticity comes at higher 
frequencies. 


The spin-down limit on strain defined in Equation 20 for known pulsars 
requires knowing the frequency fqw, its first derivative faw and the distance d 
to the star. There are other neutron stars for which no pulsations are observed, 
hence for which neither fa@w nor faw is known, but for which the distance 
and the age of the star are known with some precision. For such stars one can 
define an “age-based” limit — under the assumption of gravitar behavior since 
the neutron star’s birth in a supernova event. Using Equation 7 and a braking 
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Fig. 8 Maximum allowed (optimistic) all-sky-search ranges for gravitars vs GW frequencies 
for different assumed spin-down derivatives for Advanced LIGO sensitivity (solid magenta 
curves) and corresponding ranges for ETC sensitivity (dashed green curves). Known pulsar 
distances are shown vs the expected GW frequencies, where red dots indicate pulsars with 
accessible spin-down limits for Advanced LIGO sensitivity, and smaller green dots indicated 
pulsars with accessible spin-down limits for ETC sensitivity. These all-sky search ranges 
assume an optimistic sensitivity depth of 50 Hz—1/? (see section 3.8). 


index of 5 for mass quadrupole radiation gives the gravitar age: 


= f rot 

4 f rot 
Therefore, if one knows the distance and the age of the star, e.g., from the 
expansion rate of its visible nebula, then under the assumption that the star 
has been losing rotational energy since birth primarily due to gravitational 


wave emission, then one has the following frequency-independent age-based 
limit on strain (Wette et al., 2008): 


Rage = (2.3 x 10774) (—=") (— *) (=) (28) 
r Tt i) 


(27) 


Tgravitar = 


26 Keith Riles 


along with a corresponding frequency-dependent but distance-independent el- 
lipticity upper limit (Wette et al., 2008): 


2 
wma tannin (SEF) Y(MP)(2) oe 


The corresponding calculation for r-mode emission leads to the age-based 
strain limit relation (Owen, 2010): 


= _ 1kpe 1000 yr L, 
r—mode __ 24 ZL 
mgr aacio™ (2282) [CY (Fa), an 


along with a corresponding frequency-dependent but distance-independent r- 
mode amplitude upper limit (Wette et al., 2008): 


1000 yr\ /? (100 Hz\? 
ange = 0.076 ( ox) (=) (31) 


Yet another empirically determined strain limit upper limit can be defined 
for accreting neutron stars in binary systems, such as Scorpius X-1. The X- 
ray luminosity from the accretion is a measure of mass accumulation rate at 
the surface. As the material rains down on the surface it can add angular 
momentum to the star, which in equilibrium may be radiated away in gravi- 
tational waves. Hence one can derive a torque-balance limit (Wagoner, 1984; 
Papaloizou and Pringle, 1978; Bildsten, 1998) in the form (Watts et al., 2008): 


3 1 
R tf Mo \# 
Rtorque © (3 x 107-77) (Ex) (42) 


at 2 
1 Hz \ ? 
- ( 000 *) Fy , . (32) 
Frot 10-8 erg/cm*/s 


where F;. is the observed energy flux at the Earth of X-rays from accretion, 
M is the neutron star mass and R its radius. Taking nominal values of R 
= 10 km, M = 1.4Mo and reformulating in terms of the gravitational wave 
frequency faw (benchmarked to 600 Hz), one obtains: 


1 3 
600 Hz \ ? F 
—27 x 
ie ee Ae ( faw ) (= 2) = 


Equations 32-33 assumes the radius at which the accretion torque is ap- 
plied is the stellar surface. If one assumes the torque lever arm is the Alfvén 
radius because of the coupling between the stellar rotation and the magneto- 
sphere, then the implied equilibrim strain is ~2.4 times higher (Abbott et al., 
2019e). This limit is independent of the distance to the star. In general, varia- 
tions in accretion inferred from X-ray flux fluctuations suggest similar (slower) 
fluctuations in the equilibrium frequency, which could degrade GW detection 
sensitivity for coherent searches that assume exact equilibrium. A first at- 
tempt to address these potential frequency fluctuations for Scorpius X-1 may 
be found in (Mukherjee et al., 2018). 
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2.1.8 Assessing potential sources of neutron star non-axisymmetry 


From the above, it is clearly possible for neutron stars in our galaxy to produce 
continuous gravitational waves detectable by current ground-based detectors, 
but is it lékely that putative emission mechanisms are strong enough to give 
us a detection in the next few years. Let’s look more critically at those mech- 
anisms °. 

Isolated neutron stars may exhibit intrinsic non-axisymmetry from resid- 
ual crustal deformation (e.g., from “starquakes” due to cooling & cracking 
of the crust (Pandharipande et al., 1976; Kerin and Melatos, 2022) or due 
to changing centrifugal stress induced by stellar spin-down (Ruderman, 1969; 
Baym et al., 1969; Fattoyev et al., 2018; Giliberti and Cambiotti, 2022)), from 
non-axisymmetric distribution of magnetic field energy trapped beneath the 
crust (Zimmermann, 1978; Cutler, 2002) or from a pinned neutron super- 
fluid component in the star’s interior (Jones, 2010; Melatos et al., 2015). 
See (Haskell et al., 2015; Singh et al., 2020) for a discussion of emission 
from magnetic and thermal “mountains” and (Lasky, 2015; Glampedakis and 
Gualtieri, 2018) for recent, comprehensive reviews of GW emission mechanisms 
from neutron stars). 

Maximum allowed asymmetries depend on the neutron star equation of 
state (Johnson-McDaniel and Owen, 2013; Krastev et al., 2008) and on the 
breaking strain of the crust. Detailed molecular dynamics simulations bor- 
rowed from condensed matter theory have suggested in recent years that the 
breaking strain may be an order of magnitude higher than previously thought 
feasible (Horowitz and Kadau, 2009; Caplan and Horowitz, 2017). Analytic 
treatments (Baiko and Chugunov, 2018) indicate, however, that anisotropy 
may be important and caution that simulations based on relatively small num- 
bers of nuclei may not capture effects due to a polycrystalline structure in the 
crust. A recent cellular automaton-based simulation (i<erin and Melatos, 2022) 
of a spinning-down neutron star used nearest-neighbour tectonic interactions 
involving strain redistribution and thermal dissipation. That study found the 
resulting annealing led to emitted gravitational strain amplitudes too low to 
be detected by present-generation detectors. 

A recent revisiting of the mountain-building scenario (Gittins et al., 2020) 
finds systematically lower ellipticities to be realistic. It is argued in (Woan 
et al., 2018) that a possible minimum ellipticity in millisecond pulsars may 
arise from asymmetries of buried internal magnetic field B; (Cutler, 2002; 
Lander et al., 2011; Lander, 2014) of order of (Woan et al., 2018) 


<B;> <H.> 
~ 107° : = 34 
° (S23) (Sez). eo 
where H, is the lower critical field for superconductivity (protons in the stellar 
core are assumed to form a Type II superconductor). Hence, a buried toroidal 


6 This section draws, in part, from a recent review (Glampedakis and Gualtieri, 2018) of 
gravitational wave emission from neutron stars. 


28 Keith Riles 


(equatorial) field of ~10!1 G could yield an ellipticity at the 10~° level. It 
has been argued, on the other hand, that an explicit model of braking dy- 
namics with non-axisymmetry due to magnetic field non-axisymmetry leads 
to still smaller ellipticities, based on observed braking indices of younger pul- 
sars (de Araujo et al., 2016, 2017), where the magnetic contribution to the 
ellipticity depends quadratically on the field strength (Bonazzola and Gour- 
goulhon, 1996; Konno et al., 1999; Regimbau and de Freitas Pacheco, 2006). 
An analysis (Osborne and Jones, 2020) of internal magnetic field contributions 
to non-axisymmetric temperature distributions finds that high field strengths 
(> 10'3 G) are needed in an accreting system for GW emission to halt spin- 
up from the accretion, four orders of magnitude higher than is expected for 
surface fields in LMXBs. 

r-modes (mass current quadrupole, see section 2.1.2) offer an intriguing 
alternative GW emission source (Mytidis et al., 2015). Serious concerns have 
been raised, (Arras et al., 2003; Glampedakis and Gualtieri, 2018) however, 
about the detectability of the emitted radiation for young isolated neutron 
stars, for which mode saturation appears to occur at low r-mode amplitudes 
because of various dissipative effects (Owen, 2010). Another study, (Alford and 
Schwenzer, 2014) though, is more optimistic about newborn neutron stars. The 
same authors, on the other hand, find that r-mode emission from millisecond 
pulsars is likely to be undetectable by Advanced LIGO (Alford and Schwenzer, 
2015). 

The notion of a runaway rotational instability was first appreciated for 
high-frequency f-modes, (Chandrasekhar, 1970; Friedman and Schutz, 1978) 
(Chandrasekhar-Friedman-Schutz instability), but realistic viscosity effects seem 
likely to suppress the effect in conventional neutron star production (Lindblom 
and Detweiler, 1977; Lindblom and Mendell, 1995). Moreover, (Ho et al., 2019) 
set limits on the r-modes amplitude a for J0952—0607 below 10° based on 
the absence of heating observed in its X-ray spectrum, despite its high ro- 
tation frequency (707 Hz) which places it in the nominal r-modes instability 
window. Similarly, (Boztepe et al., 2020) set limits on a as low as 3 x 107°, 
based on observations of two other millisecond pulsars (PSR J1810+1744 and 
PSR J2241—5236) which also sit in the instability window. Another potential 
source of r-modes dissipation is from the interaction of “ordinary” and su- 
perfluid modes, leading to a stabilization window for LMXB stars (Gusakov 
et al., 2014; Kantor et al., 2020). The f-mode stability could play an important 
role, however, for a supramassive or hypermassive neutron star formed as the 
remnant of a binary neutron star merger (Doneva et al., 2015). 

In addition, as discussed below, a binary neutron star may experience direct 
non-axisymmetry from non-isotropic accretion (Owen, 2005; Ushomirsky et al., 
2000; Melatos and Payne, 2005) (also possible for an isolated young neutron 
star that has experienced fallback accretion shortly after birth), or may exhibit 
r-modes induced by accretion spin-up. 

Given the various potential mechanisms for generating continuous gravi- 
tational waves from a spinning neutron star, detection of the waves should 
yield valuable information on neutron star structure and on the equation of 
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state of nuclear matter at extreme pressures, especially when combined with 
electromagnetic observations of the same star. 

The notion of gravitational wave torque equilibrium is potentially impor- 
tant, given that the maximum observed rotation frequency of neutron stars 
in LMXBs is substantially lower than one might expect from calculations of 
neutron star breakup rotation speeds (~1400 Hz) (Cook et al., 1994). It has 
been suggested (Chakrabarty et al., 2003) that there is a “speed limit” due to 
by gravitational wave emission that governs the maximum rotation rate of an 
accreting star. In principle, the distribution of frequencies could have a quite 
sharp upper frequency cutoff, since the angular momentum emission is pro- 
portional to the 5th power of the frequency for mass quadrupole radiation. For 
example, for an equilibrium frequency corresponding to a particular accretion 
rate, doubling the accretion rate would increase the equilibrium frequency by 
only about 15%. For r-mode GW emission, with a braking index of 7, the 
cutoff would be still sharper. 

Note, however, that a non-GW speed limit may well arise from interaction 
between the neutron star’s magnetosphere and an accretion disk (Ghosh and 
Lamb, 1979; Haskell and Patruno, 2011; Patruno et al., 2012). It has also been 
argued (Ertan and Alpar, 2021) that correlation between the accretion rate 
and the frozen surface dipole magnetic field resulting from Ohmic diffusion 
through the neutron star crust in the initial stages of accretion in low mass 
X-ray binaries can explain a minimum rotation period well above the naive 
expectation. 

A number of mechanisms have been proposed by which the accretion leads 
to gravitational wave emission in binary systems. The simplest is localized 
accumulation of matter, e.g., at the magnetic poles (assumed offset from the 
rotation axis), leading to a non-axisymmetry. One must remember, however, 
that matter can and will diffuse into the crust under the star’s enormous 
gravitational field. This diffusion of charged matter can be slowed by the also- 
enormous magnetic fields in the crust, but detailed calculations (Vigelius and 
Melatos, 2010) indicate the slowing is not dramatic. Relaxation via thermal 
conduction is considered in (Suvorov and Melatos, 2019). 

Another proposed mechanism is excitation of r-modes in the fluid interior of 
the star, (Andersson, 1998; Bildsten, 1998; Friedman and Morsink, 1998; Owen 
et al., 1998) with both steady-state emission and cyclic spin-up/spin-down pos- 
sible (Levin, 1999; Heyl, 2002; Arras et al., 2003). Intriguing, sharp lines con- 
sistent with expected r-mode frequencies were reported in the accreting mil- 
lisecond X-ray pulsar XTE J1751—305 (Strohmayer and Mahmoodifar, 2014a) 
and in a thermonuclear burst of neutron star 4U 1636—536 (Strohmayer and 
Mahmoodifar, 2014b). The inconsistency of the observed stellar spin-downs for 
these sources with ordinary r-mode emission, however, suggests that a different 
type of oscillation is being observed (Andersson et al., 2014) or that the puta- 
tive r-modes are restricted to the neutron star crust and hence gravitationally 
much weaker than core r-modes (Lee, 2014). Another recent study (Patruno 
et al., 2017) suggests that spin frequencies observed in accreting LMXB’s are 
consistent with two sub-populations, where the narrow higher-frequency com- 
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ponent (~575 Hz with standard deviation of ~30 Hz) may signal an equilib- 
rium driven by gravitational wave emission. It has been suggested (Haskell 
and Patruno, 2017) that the transitional millisecond pulsar PSR J1023+0038 
(for which spin-down has been measured in both accreting and non-accreting 
states) shows evidence for mountain building (or r-modes) during the accretion 
state, based on different spin-downs observed in accreting vs. non-accreting 
states. It has also been argued (Bhattacharyya, 2020) that J1023+0038 shows 
evidence for a permanent ellipticity in the range 0.48 — 0.93 x 10~°. An anal- 
ysis (Chen, 2020) of three transitional millisecond pulsars and ten redbacks 
concluded their ellipticities ranged over 0.9 — 23.4 x 1079. 

A recent analysis (De Lillo et al., 2022) based on the absence of evidence of 
a stochastic gravitational wave background emitted by a population of neutron 
stars with a rotational frequency distribution similar to that of known pulsars 
inferred that the average ellipticity of the galactic population is less than 
~2 x 1078. 


2.1.4 Particular GW targets 


In the following, particular neutron star targets for gravitational wave searches 
are discussed in the following categories: known young pulsars with high spin- 
down rates; known high-frequency millisecond pulsars; neutron stars in super- 
nova remnants, neutron stars in low-mass X-ray binary systems; and particular 
directions on the sky. 


2.1.4.1 Known young pulsars with high spin-down rates 


A young pulsar with a high spin-down rate presents an attractive target. Its 
age offers the hope of a star not yet annealed into smooth axisymmetry, a hope 
strengthened by the prevalence of observed timing glitches among young stars. 
A high spin-down rate not only makes it more likely that the spin-down limit 
is accessible, but also suggests a star with a reservoir of magnetic energy, some 
of which could give rise to non-axisymmetry. From the Advanced LIGO O1, 
O2 and O3 data sets more than 20 pulsars were spin-down accessible (Abbott 
et al., 2019b; Abbott et al., 20211) (see section 4.1), but most correspond to 
ellipticities of ~(10~4-10~3). A small number are highlighted here, for which 
ellipticities below 10~° are accessible already or with a 2-year data run at 
Advanced LIGO design sensitivity (“O4/05 run”). 


— Crab (PSR J0534+2200) — This pulsar, created in a 1054 A.D. super- 
nova observed by Chinese astronomers and discovered in 1968 (Staclin and 
Reifenstein, 1968), has received more attention from LIGO / Virgo analysts 
than any other. Its spin-down limit was first beaten in the initial LIGO 
data set S55 (Abbott et al., 2008b), and now has been beaten (O03 data) by 
a factor of ~100 (Abbott et al., 20211) (see section 4.1), leading to a 95% 
upper limit on ellipticity of 1.0 x 10~°. For a 2-year 04/05 run, this sensi- 
tivity reaches ~2 x 10~°. Spinning at just below 30 Hz, its nominal faw is 
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just below 60 Hz, making the detector spectrum susceptible to power mains 
contamination (including non-linear upconversion, see section 3.7) in the 
LIGO interferometers. Its inferred rotational kinetic energy loss rate based 
on its spin-down is dE/dt ~ —5 x 10°8 erg s~+, assuming I, = 10°° kg m? 
(10° ¢ cm?). 

— Vela (PSR J0835—4510) — Although older and lower in frequency than 
the Crab with a higher ellipticity spin-down limit (1.9 1073), the Vela pul- 
sar, discovered in 1968 (Large ct al., 1968), is nonetheless interesting, given 
its frequent glitches (Manchester, 2018; Ashton et al., 2019). Its 04/05 el- 
lipticity sensitivity reaches ~8 x 10~®. Spinning at just above 11 Hz, its 
nominal few is about 22 Hz, where detector noise is several times higher 
than at the Crab frequency. Its inferred dE /dt ~ —7 x 10°° erg s7!. 

— PSR J0537—6910 — This pulsar, observed to pulse only in X-rays, is dis- 
tant (~50 kpc in the Large Magellanic Cloud). With a rotation frequency 
of ~62 Hz, its nominal GW frequency of 124 Hz is quite high for a young 
pulsar (magnetic dipole spin-down age ~ 5000 yr), and its spin-down en- 
ergy loss is comparable to the Crab’s. It is also extremely glitchy (~1 per 
100 days) (Antonopoulou et al., 2018; Ferdman et al., 2018) and as noted 
above, may show evidence of r-mode emission between glitches (Andersson 
et al., 2018; Ho et al., 2020) (which would imply a GW frequency at ~90 
Hz. Its inferred dE /dt ~ —5 x 108 erg s7?. 

— PSR J1400—6325 — This relatively recently discovered X-ray pulsar (Re- 
naud et al., 2010) lies in a supernova remnant 7-10 kpc away and displays 
a spin-down energy about 1/10 of the Crab pulsar’s, but may be younger 
than 1000 yr. With a spin frequency of ~32 Hz, its nominal few is 64 Hz, 
comparable to the Crab’s, but farther from the troublesome 60-hz power 
mains. Its inferred dE/dt ~ —5 x 10°” erg s~. 

— PSR J1813—1749 — First detected as a TeV y-ray source (Aharonian 
et al., 2005), this star was found to exhibit non-thermal X-ray emission 
and to have a tentative association with a radio supernova remnant G12.8- 
0.0 (Brogan et al., 2005) suggesting a distance greater than 4 kpc and an 
age perhaps younger than 1000 years. X-ray pulsations detected still later 
with a period of 44 ms confirmed a pulsar source and posited an association 
with a young star cluster at 4.7 kpc (Gotthelf and Halpern, 2009)], while 
yielding a nominal pulsar spin-down age of 3.3-7.5 kyr. A more recent 
detection of highly dispersed radio pulsations, however, suggest a distance 
of 6 or 12 kpc (Camilo et al., 2021), depending on electron dispersion model, 
casting doubt on the association with the star cluster. The spin frequency 
of 22 Hz yields a nominal GW frequency of ~45 Hz, and the frequency 
derivative imply dE/dt ~ —6 x 10°” erg s~. 


2.1.4.2 Known high-frequency millisecond pulsars 
Because nearly all millisecond pulsars are old, with some characteristic ages 


greater than 10 billion years, they can be assumed to retain little asymmetry 
from their initial formation or from the accretion that spun them up. Thus 
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one sees low spin-down for this population in Fig. 5 and hence low inferred 
maximum ellipticities in Fig. 6. On the other hand, the vast energy reser- 
voirs in their rotation and the quadratic dependence of hp on frequency still 
makes these stars potentially intriguing. As noted above, there may be em- 
pirical evidence for a minimum ellipticity of order ~10~® (Woan et al., 2018). 
Highlighted below are particular millisecond pulsars of interest in the coming 
years. 


— PSR J0711—6830 This isolated star at a distance of 0.11 kpc, with a 
nominal faw ~ 364 Hz, a spin-down upper limit of 1.2 x 10~?° and cor- 
responding maximum ellipticity of 9.4 x 10~°, is the first MSP to have its 
spin-down limit beaten (in early O3 data, see section 4.1). 

— PSR J0437—4715 — This binary star at a distance of 0.16 kpc, with 
a nominal fgw ~ 347 Hz, a spin-down upper limit of 7.8 x 10727 and 
corresponding maximum ellipticity of 9.7 x 10~° also had its spin-down 
limit beaten (in the full O3 data). 

— PSR J1737—0811 This binary star at a distance of 0.21 kpc, with a 
nominal faw ~ 479 Hz, a spin-down upper limit of 5.3 x 1072? and corre- 
sponding maximum ellipticity of 4.6 x 10~°, will likely have its spin-down 
limit beaten by the 04/05 data set. 

— PSR J1231—1411 This binary star at a distance of 0.42 kpc, with a 
nominal faw ~ 543 Hz, a spin-down upper limit of 2.8 x 1072? and corre- 
sponding maximum ellipticity of 3.8 x 10~°, will likely have its spin-down 
limit beaten by the 04/05 data set. 

— PSR J2124—3358 This binary star at a distance of ~0.4 kpc, with a 
nominal faw ~ 406 Hz, a spin-down upper limit of 2.3 x 10727 and corre- 
sponding maximum ellipticity of 5.6 x 10~°, will likely have its spin-down 
limit beaten by the 04/05 data set. 

— PSR J1643—1224 This binary star at a distance of 0.79 kpc’, with a 
nominal fgw ~ 433 Hz, a spin-down upper limit of 2.1 x 107?? and cor- 
responding maximum ellipticity of 8.0 x 10-°, may not have its spin-down 
limit beaten by the 04/05 data set, but as noted by (Woan et al., 2018), 
would have the highest GW SNR of any known star if its ellipticity were 
107°. 


2.1.4.8 Central compact objects and Fomalhaut b 


Not every neutron star of interest has been detected to pulsate. Central 
compact objects (CCOs) at the heart of supernova remnants present especially 
intriguing targets, especially those in remnants inferred from their size and ex- 
pansion rate to be young (De Luca, 2008). There may be direct evidence of a 
neutron star, such as from thermal X-rays emitted from a hot surface or from 
X-rays due to interstellar accretion, or there may be indirect evidence from 
a pulsar wind nebula driven by a fast-spinning star at the core. Most GW 


7 A recent parallax measurement (Reardon et al., 2021), though, finds a distance for 
J1643—1224 of 1.2103 kpc. 
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searches to date for a CCO lacking detected pulsations have focused on the 
particularly promising source, Cassiopeia A, but in recent years, such searches 
have also been carried out for as many as 15 supernova remnants (Abbott et al., 
2019a; Abbott et al., 2021k). Highlighted below are particular supernova rem- 
nants (“G” naming terminology based on the Green Catalog (Green, 2014), 
see also (Ferrand and Safi-Harb, 2012)) with known or suspected central com- 
pact objects, in addition to an object, Fomalhaut b, originally thought to be 
an exoplanet, but which may be a nearby neutron star. Results from searches 
for these targets are presented further below in section 4.2. 


— Cassiopeia A — Cas A (G111.7—2.1) is perhaps the most promising exam- 
ple of gravitational wave CCO source in a supernova remnant. Its birth af- 
termath may have been observed by Flamsteed (Hughes, 1980) ~340 years 
ago in 1680, and the expansion of the visible shell is consistent with that 
date (Fesen ct al., 2006). Hence Cas A, which is visible in X-rays (Tanan- 
baum, 1999; Ho et al., 2021) but shows no pulsations (Halpern and Got- 
thelf, 2009). is almost certainly a very young neutron star at a distance of 
about 3.3 kpc (Reed et al., 1995; Alarie et al., 2014). From equation 28, one 
finds an age-based strain limit of ~1.2 x 10724, which is readily accessible 
to LIGO and Virgo detectors in their most sensitive band. 

— Vela Jr. — This star (G266.2—1.2) is observed in X-rays (Pavlov et al., 
2001; Kargaltsev et al., 2002; Becker et al., 2006) and is potentially quite 
close (~0.2 kpc) and young (690 yr) (Iyudin et al., 1998), but searches 
have also conservatively assumed more a more pessimistic distance (0.9 
kpc) and age (5100 yr), based on other measurements (Allen et al., 2015). 
The optimistic age and distance assumptions lead to an age-based strain 
limit of ~1.4 x 107%, even more accessible than the Cas A limit. Even the 
pessimistic age-base limit of 1.1 x 1074 is only slightly lower than that of 
Cas A. It has been argued (Ming ct al., 2016a) that a search over multiple 
CCOs, optimized for most likely detection success given fixed computing 
resources, favors focusing those resources on Vela Jr. over other CCOs, 
including Cas A. 

— G347.3—0.5 — An X-ray source (Slane et al., 1999; Ho et al., 2021) is 
consistent with the core of this supernova remnant, the nearness (~0.9 
kpc) and youth (1600 yr) of which make a search aimed at the remnant’s 
center intriguing, as they yield an age-based strain limit of ~2.0 x 10724 — 
higher than that of Cas A. 

— G1.9+0.3 This supernova remnant, the youngest in the galaxy at 100 
yr (Reynolds et al., 2008), has no detected CCO at its core, which is con- 
sistent with a Type IA supernova’s having left no neutron star behind. 
Nonetheless, its youth make it interesting despite this doubt and its dis- 
tance (8.5 kpc), yielding an age-based strain limit of ~8.4 x 107°. 

— Fomalhaut b This object was assumed to be an extrasolar planet (Kalas 
et al., 2008) until (Neuhduser et al., 2015) noted that the absence of de- 
tected infrared radiation could indicate the object is a remarkably nearby 
neutron star (~0.01 kpc). The absence of attempted X-ray detection with 
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Chandra observations (Poppenhaeger et al., 2017), however, disfavors its 
being a young, hot neutron star. More recent evidence (Gaspar and Ricke, 
2020) argues, in fact, that the optical observations point to a planetesimal 
collision). 


2.1.4.4 Neutron stars in low-mass X-ray binary systems 


Because of its high X-ray flux (F, ~ 3.9 x 1077 erg cm7? s~! (Watts 
et al., 2008)) and the torque-balance relation for low-mass X-ray binaries [equa- 
tion (32)], Scorpius X-1 is thought to be an especially promising search target 
for advanced detectors and has been the subject of multiple searches in ini- 
tial and Advanced gravitational wave detector data. From equation (33), one 
expects a strain amplitude limited by (Abbott et al., 2007a; Messenger et al., 
2015) 


600 Hz \ 2 
a) , (35) 


how (3x 107) ( 7 
GW 


While Sco X-1’s rotation frequency remains unknown (Galaudage et al., 2021), 
its orbital period is well measured, (Gottlieb et al., 1975; Wang et al., 2018) 
which allows substantial reduction in search space. A similar but less bright 
LMXB system is Cygnus X-2 (Premachandra et al., 2016) at a distance of 7 
kpc and an average flux F, = 11 x 10~® erg/cm? s~+ (Galloway et al., 2008), 
yielding a torque-balance strain limit about 20 times lower than that of Sco 
X-1. Unlike Sco X-1 which is assumed but not known to contain a neutron star 
(as opposed to a black hole with an accretion disk), Cyg X-2 has displayed 
thermonuclear bursts, confirming the presence of a neutron surface. 

Another interesting class contains “accreting X-ray millisecond pulsars” 
(AXMPs) which are fast-spinning neutron stars in LMXBs that show spo- 
radic outbursts during accretion episodes (when “active”) from which rotation 
frequencies can be determined. When active, the frequencies can increase or 
decrease, while frequencies between outbursts (when “quiescent”) generally 
decrease. One could hope to detect CW radiation from either active or qui- 
escent phases. Although the limited durations of bursts and their stochastic 
nature constrain potential search sensitivity, it is during such outbursts when 
one might expect the largest generation of non-axisymmetries or excitation 
of r-modes. The fastest-spinning stars, such as IGR J00291+5934 at fro ~ 
599 Hz and a distance of ~4 kpc (Torres et al., 2008; Patruno, 2017), offer 
deeper probing of equatorial ellipticity and r-mode amplitude. Current search 
sensitivities to strain amplitude (Abbott et al., 2022d) remain an order of 
magnitude or more away from inferred spin-down limits (~10~78-10~2”), but 
improvements in detector sensitivity, search methodology and potential future 
electromagnetic observations make this type of source potentially intriguing 
in the coming years. 


2.1.4.5 Particular sky directions 


Searches for Continuous-Wave Gravitational Radiation 35 


In addition to known (or suspected) neutron stars, there are other local- 
ized sky regions or points where a directed search might yield a continuous 
gravitational wave detection. Listed below are possibilities that have attracted 
attention in recent years. 


— Galactic center — The vicinity of the galactic center (Sgr A*) is particu- 
larly interesting(Aasi et al., 2013a), as an active, star-forming region with 
known pulsars (Deneva et al., 2009). Moreover, it is highly likely that only 
a small fraction of pulsars near the galactic center have been detected to 
date, since there is extreme dispersion and scattering of radio signals along 
the propagation line to the Earth (Lazio and Cordes, 1998). The infer- 
ence of there being many hidden pulsars is supported by ~20 pulsar wind 
nebula candidates detected within 20 pc of Sgr A* (Muno et al., 2008). 
In addition, searches for dark-matter annihilation signals have detected an 
excess of high-energy gamma ray emission from the galactic center region 
above what is expected from conventional models of diffuse gamma-ray 
emission and catalogs of known gamma-ray sources (Ackermann et al., 
2017), a tension which may be resolved by the existence of a hidden popu- 
lation of millisecond pulsars (Abazajian, 2011). A systematic radio survey 
of the central 1 parsec of Sgr A* at a frequency of 15 GHz (Macquart 
and Kanekar, 2015), high enough to reduce dispersion and scattering sub- 
stantially, yielded no detections, but the rapidly falling spectrum of most 
pulsars makes detection at 15 GHz at that distance difficult. This survey 
obtained a 90% CL upper limit of 90 on the number of pulsars within 1 par- 
sec of Sgr A*, assuming the population there is similar to known pulsars. 
Unfortunately, the ~8.5 kpc distance to the galactic center makes CW 
searches challenging with present detector sensitivities. Only stars with 
extreme ellipticities are accessible to Advanced LIGO at design sensitivity 
(see Figure 8). At the same time, however, young neutron stars are those 
most likely to exhibit such ellipticities. 

— Globular cluster cores — One normally associates globular clusters with 
ancient stellar populations and might expect, at best, to see only pulsars 
that are themselves ancient — recycled and well annealed millisecond pul- 
sars. Indeed many MSPs are seen in globular clusters (Freire, 2012). For 
example, Tucanae 47 is known to host at least 25 MSPs (Freire et al., 
2017). Nonetheless, not all observed pulsars in globular clusters seem to be 
old (Freire, 2012). A plausible explanation is that the dense core of a glob- 
ular cluster leads to multibody exchange interactions in which a previously 
recycled but decoupled neutron star acquires a close new companion that 
proceeds to overflow its Roche lobe, leading to new accretion. Another, 
related mechanism is possible debris accretion triggered by multibody in- 
teractions, given that some pulsars are known to host debris disks and even 
planets (Abbott et al., 2017m). The well localized core of a globular cluster 
makes a deep, directed search tractable. 

— High-latitude Fermi sources The Fermi satellite’s LAT experiment has 
detected ~100 previously unknown gamma ray pulsars since observing be- 
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gan in 2008. Gamma ray pulsars tend to be sources with low variability 
and relatively low spectral cutoffs, and most lie near the galactic plane, 
as expected. Fermi-LAT point sources well outside the galactic plane tend 
to be extagalactic, e.g., active galactic nuclei, but an intriguing possib- 
lity is that a source with high galactic latitude could be a galactic neu- 
tron star, in which case the high latitude favors a nearby source (Sanders, 
2016), consistent with a scale height of ~600 pc with respect to the galactic 
plane observed for known pulsars (Lyne and Graham-Smith, 2006). Argu- 
ing against this possibility, however, are extensive searches for gamma-ray 
pulsations from pulsar-like Fermi-LAT sources (see, e.g.., (Clark et al., 
2018)), based in part on algorithms developed for CW gravitational wave 
searches (Pletsch et al., 2011). On the other hand, such searches are chal- 
lenged to probe binary sources with large accelerations, suggesting that 
CW searches directed at such sources include algorithms sensitive to bi- 
nary sources (Neunzert, 2019). 


In between all-sky searches and directed searches for single sky points reside 
“spotlight” searches, in which a patch of sky is searched more deeply than in 
all-sky searches (with increased computational cost), but less deeply than is 
computationally feasible for a single sky location. Such spotlights have been 
applied in searches for a broad star-forming region along two directions of the 
Orion spur of the local galactic spiral arm (Aasi ect al., 2016b) and toward 
the galactic center region, including the globular cluster Terzan 5 (Dergachev 
et al., 2019). 


2.2 Axion clouds bound to black holes 


An intriguing potential connection between gravitational waves and the still- 
unknown missing dark matter of the Universe comes from the possibility that 
the dark matter is composed of ultralight, electromagnetically invisible bosons, 
such as axions. One novel idea is that these bosons could be disproportionately 
found in the vicinity of rapidly spinning black holes (Arvanitaki et al., 2010; 
Arvanitaki and Dubovsky, 2011). The ultralight particles could, in principle, 
be spontaneously created via energy extraction from the black hole’s rota- 
tion (Penrose, 1969; Christodoulou, 1970) and form a Bose-Einstein “cloud” 
with nearly all of the quanta occupying a relatively small number of energy 
levels. For a cloud bound to a black hole, the approximate inverse-square law 
attraction outside the Schwarzchild radius (rgchwarz. = aout) leads to an 
energy level spacing directly analogous to that of the hydrogen atom (Arvani- 
taki et al., 2010; Baumann et al., 2019). The number of quanta occupying the 
low-lying levels can be amplified enormously by the phenomenon of superradi- 
ance in the vicinity of a rapidly spinning black hole (with angular momentum 
that is a signficant fraction of the maximum value allowed in General Rela- 
tivity). The bosons in a non-s (¢ > 0) negative-energy state can be thought 
of as propagating in a well formed between an ¢-dependent centrifugal barrier 
at r > Tgchwarz, and a potential rising toward zero as r — oo; wave function 
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penetration into the black hole ergosphere permits transfer of energy from the 
black hole spin (Zel’Dovich, 1971; Misner, 1972; Starobinskii, 1973) into the 
creation of new quanta. 

Two particular gravitational wave emission modes of interest here can 
arise in the axion scenario, both potentially leading to intense coherent ra- 
diation (Arvanitaki et al., 2015). In one mode, axions can annihilate with each 
other to produce gravitons with frequency double that corresponding to the 
axion mass: feraviton = WMaxionC /h. In another mode, emission occurs from 
level transitions of quanta in the cloud. This Bose condensation is most pro- 
nounced when the reduced Compton wavelength of the axion is comparable 
to but larger than the scale of the black hole’s Schwarzchild radius: 


r h > 2 GMpu 
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where fi is the reduced Planck constant and G' is Newton’s gravitational con- 
stant. A key parameter governing detectability is a parameter analogous to 
the electromagnetic fine structure constant: 


GmMaxionMpu 
hie , 


a (38) 
where both the growth rate of a cloud upon black hole formation and the 
amplitude of gravitational wave emission due to axion annihilation depend 
on high powers of a. Hence small a@ impedes detection; at the same time, 
superradiance itself requires (Isi et al., 2019): 


-1 
a< mx (1+ VI-x?) < - (39) 


where y is the dimensionless black hole spin proportional to its total angular 
momentum magnitude J: x = ape. and m is the quantum number corre- 
BH 


sponding to the axion’s orbital angular momentum projection along the spin 
axis of the black hole (the first level to be populated in a newborn black hole is 
m = 1 (isi et al., 2019)). Hence the range of a (and therefore axion mass) for 
which a particular black hole produces superradiance may be narrow. In gen- 
eral, more massive black holes produce stronger signals over wider ranges in 
axion mass. Clouds composed of ultralight vector or tensor bosons would lead 
to stronger, but shorter-lived signals (Siemonsen and East, 2020; Brito et al., 
2020). Nominal limits on axion masses can be placed based on the existence 
of high-spin binary black holes in our galaxy (Arvanitaki et al., 2017; Cardoso 
et al., 2018), but those limits are subject to uncertainties in inferred black 
hole spins (Reynolds, 2014; McClintock et al., 2014) and may be invalidated 
by tidal disruption effects from the companion star (Cardoso et al., 2020). 
Constraints have also been inferred from spin measurements in the population 
of binary black hole merger detections (Ng et al., 2021). 
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Given the many orders of magnitude of uncertainty in, for example, axion 
masses that could account for dark matter(Bertone, 2010), the relatively nar- 
row mass window accessible to currently feasible CW searches (1-2 orders of 
magnitude) makes searching for such an emission a classic example of “lamp- 
post” physics, where one can only hope that nature places the axion in this 
lighted area of a vast parameter space. 

In principle, searching for these potential CW sources requires no funda- 
mental change in the search methods described below, but search optimization 
can be refined for the potentially very slow (and positive) frequency evolution 
expected during annihilation emission (as the relative magnitude of the axion 
field’s binding energy decreases). In addition, for a known black hole location, a 
directed search can achieve better sensitivity than an all-sky search. For string 
axiverse models, however, the axion cloud (Arvanitaki et al., 2010; Arvanitaki 
and Dubovsky, 2011; Yoshino and Kodama, 2014, 2015) can experience signif- 
icant self-interactions which can lead to appreciable frequency evolution of the 
signal and to uncertainty in that evolution, a complication less important for 
the postulated QCD axion (Arvanitaki ct al., 2015). In an optimistic scenario 
with many galactic black holes producing individually detectable signals, (Zhu 
et al., 2020) points out that the signals would all lie in a very narrow band, 
complicating CW searches, which typically implicitly assume no more than 
one detectable signal in narrow bands. 

Until recently, most published searches have not been tailored for a black 
hole axion cloud source, but instead existing (non-optimized) limits on neutron 
star CW emission could be reinterpreted as limits on such emission (Arvan- 
itaki et al., 2015; Dergachev and Papa, 2019; Palomba et al., 2019). More 
recently, though, searches have been carried out that exploit the narrow spin- 
up parameter space expected for such sources (Sun et al., 2020; Abbott et al., 
2022b). 

One interesting suggestion includes the possibility that a black hole formed 
from the detected merger of binary black holes or neutron stars could provide a 
natural target for follow-up CW searches (Arvanitaki et al., 2017; Ghosh et al., 
2019; Isi et al., 2019). Recent studies (Brito et al., 2017b,a; Tsukada et al., 
2019, 2021) argue that the lack of detection of a stochastic gravitational radi- 
ation background from the superposition of extragalactic black holes already 
places significant limits on axion masses relevant to CW searches. Another 
recent study (isi et al., 2019) examined in detail the prospects for detecting 
superradiance from both post-merger black hole remnants and known black 
holes in galactic X-ray binaries, as Cygnus X-1. 


3 Continuous Wave Search Methods 


Being realistic, we must acknowledge that the first discovered CW signal will 
be exceedingly weak compared to the transient signals detected to date, an 
assumption borne out by many unsuccessful CW searches to date. One must 
integrate the signal over a long duration to observe it with statistical signf- 
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icance. Those long integrations in noise that is instantaneously much higher 
in amplitude require application of assumed signal templates to the data. In 
general, the more restrictive is the model, the better is the achievable signal-to- 
noise ratio, as one can search over a smaller volume of source parameter space. 
The following sections discuss the challenges faced in searching over larger pa- 
rameter space volumes, a common classification of general search methods, 
and specific algorithms devised to meet the challenges. 


3.1 Challenges in CW signal detection and types of searches 


At first glance, it may seem puzzling that a signal due to a rapidly spin- 
ning neutron star is challenging to find. One might expect a simple discrete 
Fourier transform of the data stream to reveal a sharp spike at the nomi- 
nal frequency. There are several severe complications, however, for most CW 
searches. For concreteness, imagine that a signal is weak enough to require 
a coherent, phase-preserving 1-year integration time 7T,.y. The nominal fre- 
quency resolution from a discrete Fourier transform (DFT*) is then 1/yr ~ 30 
nHz. In order for the signal’s central frequency to remain in the same DFT 
bin (integer index into the transform result, see Eqn. 52 below) during that 
year, its first derivative f would need to satisfy fTon <1/Teon, or f < 107! 
Hz/s and its second derivative f < 6 x 10~?8 Hz/s?. In practice, not only are 
Doppler modulations of detected frequency due to the Earth’s motion much 
larger than these values, as discussed below, but the frequency derivative of a 
detectable source is typically also much larger, in order for its rotational ki- 
netic energy loss to be compatible with detection (detectable spin-down limit). 
If the precise frequency evolution of the source is known already from radio 
or gamma-ray pulsar timing (assuming a fixed EM/GW phase relation), then 
one can make corrections for that evolution via barycentering, discussed below, 
without SNR degradation as long as the uncertainties in frequency derivatives 
are well below the above constraints. 

For sources with large frequency uncertainties, however, especially those 
with unknown frequencies, correcting for intrinsic source frequency evolution 
and for modulations due to the Earth’s motion incurs a substantial computing 
cost for searching over parameter space. Because of these costs, it is useful 
to categorize CW searches broadly into three categories (Prix, 2009) (while 
recognizing there are special cases that fall near the boundaries). 


1. Targeted searches in which the star’s position and rotation frequency are 
known, 7.e., known radio, X-ray or y-ray pulsars; 

2. Directed searches in which the star’s position is known, but rotation fre- 
quency is unknown, e.g., a non-pulsating X-ray source at the center of a 
supernova remnant; and 


8 CW search literature frequently refers to SFTs, “Short” discrete Fourier transforms, 
where short is relative to the span of an observing run, but which may correspond to 
coherence times as long as hours. 
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3. All-sky searches for unknown neutron stars. 


The volume of parameter space over which to search increases in large steps 
as one progresses through these categories. In each category a star can be iso- 
lated or binary. For 2) and 3) any unknown binary orbital parameters further 
increase the search volume, making a subclassification helpful, as discussed 
below. In general, the greater the a priori knowledge of sources parameters, 
the more computationally feasible it is to integrate data coherently for longer 
time periods in order to improve strain amplitude sensitivity. 

To illustrate, consider a directed search for a source of known location but 
with unknown frequency and unknown frequency derivatives, where the signal 
phase is expanded in truncated Taylor form in the source frame time 7 with 
respect to a reference time 70: 


P(r) & @94+27]|fs(7 —7)4 f(T T7)° 4 fs(v To)" |, (40) 


Using the phase evolution model of equation 40, if we wish to preserve 
phase fidelity to a tolerance A® over a coherence time Tyoh © T — To, then we 
need (in a naive estimate) to know the frequency and derivatives to a tolerance 
better than 
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Hence the numbers of steps to take in few, few, and fgw to cover a given 
range in the parameters are proportional to Tron, 2, and T3),, respectively — 
if it’s necessary to step at all in those derivatives. Naively, for a search over a 
long enough coherence time to require multiple steps in fow; one has a tem- 
plate count proportional to TS. and, presumably, pays a price proportional to 
another factor of To, in computational cost in processing the associated data 
volume. In principle, then, the computational cost of a coherent search scales 
as the 7th power of the coherence time used, although, in practice the scaling 
tends not to be as extreme because the numbers of steps needed for faw can 
be small integers that take on new discrete values only slowly with increased 
Teon- In practice, these considerations for a 2nd frequency derivative come into 
play for only directed searches or for the deep follow-up of outliers from all-sky 
searches, when segment coherence times exceed several days. Section 3.6 will 
discuss more quantitatively the placement of search templates in parameter 
space to maintain acceptable phase tolerance. 

In carrying out all-sky searches for unknown neutron stars, the computa- 
tional considerations grow worse. The corrections for Doppler modulations and 
antenna pattern modulation due to the Earth’s motion must be included, as 
for the targeted and directed searches, but the corrections are sky-dependent, 
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and the spacing of the demodulation templates is dependent upon the inverse 
of the coherence time of the search. Specifically, for a coherence time Tyon the 
required angular resolution is (Abbott et al., 2008a) 


0.5cd0f 
f [usin(9)}max’ 


where 6 is the angle between the detector’s velocity relative to a nominal source 
direction, where the maximum relative frequency shift [vsin()|max/¢ © 1074, 
and where df is the size of the frequency bins in the search. For 6f = 1/Tton, 


one obtains: 
inut H 
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where Teoh = 30 minutes has been used in several all-sky searches to date. Be- 
cause the number of required distinct points on the sky scales like 1/(60)?, the 
number of search templates scales like (Tyon)?(fs)? for a fixed signal frequency 
fs. Now consider attempting a search with a coherence time of 1 year for a 
signal frequency f, = 1 kHz. One obtains 60 ~ 0.3 wrad and a total number 
of sky points to search of ~(10'*) — again, for a fixed frequency. Adding in the 
degrees of freedom to search over ranges in f,, f, and f, (and higher-order 
derivatives, as needed) makes a brute-force, fully coherent 1-year all-sky search 
hopelessly impractical, given the Earth’s present total computing capacity. 

As a result, tradeoffs in sensitivity must be made to achieve tractability 
in all-sky searches. The simplest tradeoff is to reduce the observation time 
to a computationally acceptable coherence time. It can be more attractive, 
however, to reduce the coherence time still further to the point where the to- 
tal observation time is divided into N = Tops/Tcon, segments, each of which 
is analyzed coherently and the results added incoherently to form a detec- 
tion statistic. One sacrifices intrinsic sensitivity per segment in the hope of 
compensating (partially) with the increased statistics from being able to use 
more total data. In practice, for realistic data observation spans (weeks or 
longer), the semi-coherent approach gives better sensitivity for fixed compu- 
tational cost and hence has been used extensively in both all-sky and directed 
searches (Prix and Shaltev, 2012). One finds a strain sensitivity (threshold for 
detection) that scales approximately as the inverse fourth root of N (Abbott 
et al., 2005a). Hence, for a fixed observation time, the strain sensitivity de- 
grades roughly as N4 as Teoh decreases (see (Wette, 2012) for a discussion 
of variations from this scaling). This degradation is a price one pays for not 
preserving phase coherence over the full observation time, in order to make 
the search computationally tractable. An important virtue of semi-coherent 
searches methods, however, is robustness with respect to deviations of a signal 
from an assumed coherent model. 

In general, fully coherent search methods are potentially the most sensi- 
tive, but their applicability depends on several considerations, perhaps the 
most important being sheer computational tractability. Even when tractable 
for a particular search, moreover, a fully coherent initial search stage may incur 
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a statistical trials factor large enough to make a putative detection question- 
able, because of the necessarily ultra-fine search needed to probe coherently 
a multi-dimensional signal parameter space. Applying a priori constraints in- 
stead, when available from electromagnetic observations or theoretical expec- 
tation, can reduce the parameter space volume and hence trials factor, making 
a detection more convincing. For example, a “5 o” detection of a signal from 
a known pulsar in a targeted search might not qualify as even a weak outlier 
in an all-sky search, much less as a discovery. (Dergacheyv, 2010b) discusses 
the tradeoff between fully templated and “loose coherence” methods (see sec- 
tion 3.3.4) in a broad parameter space, arguing against brute-force template 
application. 

In the next sections, a variety of general approaches and specific algorithms 
will be presented, methods that attempt to achieve tradeoffs best suited to 
particular CW search types. 


3.2 Signal model 


CW searches must account for large phase modulations (or, equivalently, fre- 
quency modulations) of the source signal due to detector motion and poten- 
tially due to source motion (expecially for binary sources). The precision of the 
applied modulation corrections must be high in the case of targeted searches, 
which use measured ephemerides from radio, optical, X-ray or y-ray observa- 
tions valid over the gravitational wave observation time. The precision must 
also be high in following up outliers from directed or all-sky searches, while 
much less precision is needed in the first stage of hierarchical searches. This 
section describes the intrinsic signal model assumed, along with the expected 
modulations due to detector motion. 

For the Earth’s motion, one has a daily relative frequency modulation of 
Urot /C 10-® and a much larger annual relative frequency modulation of 
Vorb/C 10-4. The pulsar astronomy community has developed a powerful 
and mature software infrastructure for measuring ephemerides and applying 
them in measurements, using the TEMPO 2 program (Hobbs et al., 2006). The 
same physical corrections for the Sun’s, Earth’s and Moon’s motions (and for 
the motion of other planets), along with general relativistic effects including 
gravitational redshift in the Sun’s potential and Shapiro delay for waves pass- 
ing near the Sun, have been incorporated into the LIGO and Virgo software 
libraries (LIGO Scientific Collaboration, 2018; Astone et al., 2002a). 

Consider an isolated, rotating rigid triaxial ellipsoid (conventional model 
for a GW-emitting neutron star), for which the strain waveform detected by 
an interferometer can be written as 


1 + cos?(v) 


h(t) = Fy(t,b) ho 


cos(@(t)) + Fy (t,wW) ho cos(z) sin(®(#)), 
(46) 

where 1 is the angle between the star’s spin direction and the propagation 

direction k of the waves (pointing toward the Earth). Fy and Fy are the 
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(real) detector antenna pattern response factors (—1 < F,, Fy <1) to the + 
and x polarizations. Fy and F,, depend on the orientation of the detector and 
the source, and on the polarization angle % (Abbott et al., 2004). Here, &(t) 
is the phase of the signal, which can often usefully be Taylor-expanded as in 
equation 40, in the solar system barycenter time 7 with apparent frequency 
derivatives with respect to detector-frame time arising from source motion. 
A more general signal model with GW emission at both once and twice the 
rotation frequency is considered in (Jaranowski ct al., 1998), with effects of free 
precession addressed in (Jones and Andersson, 2001, 2002; Van Den Broeck, 
2005; Gao et al., 2020), and a convenient reparametrization is presented in 
(Jones, 2015). 

Explicitly, the time of arrival of a signal at the solar system barycenter 
(SSB), 7(t), can be written in terms of the signal time of arrival ¢t at the 
detector: : 

Yq: k 
T(t) = t+ot = ¢t : +t Ano + Aso. (47) 
where rq is the position of the detector with respect to the SSB, and Age and 
Aso are solar system Einstein and Shapiro time delays, respectively (Taylor, 
1992; Hobbs et al., 2006). 

Equation 47 implicitly assumes planar gravitational wavefronts and ne- 
glects proper motion of the source (transverse to the line of sight), corrections 
for which are common in radio pulsar astronomy (Lorimer and Kramer, 2005; 
Lyne and Graham-Smith, 2006). In principle, long-duration (multi-year) fully 
coherent observations of a near-enough (~100 pc), high-frequency (~1 kHz) 
CW source would allow inference of its distance from determination of the 
wavefront curvature (Seto, 2005). Similarly, multi-year coherent observations 
of a high-frequency source would need to account for significant proper mo- 
tions (~50 mas/year, typical of known pulsars) (Covas, 2020). Both wavefront 
curvature and proper motion have been neglected in CW searches for unknown 
sources to date because the coherence times used in the searches don’t require 
those corrections, but in the happy event of a future detection and subsequent 
extended observations, these corrections may become relevant. 

Existing gravitational wave detectors are far from isotropic in their re- 
sponse functions. In the long-wavelength limit, Michelson interferometers have 
an antenna pattern sensitivity with polarization-dependent maxima normal to 
their planes and nodes along the bisectors of the arms. As the Earth rotates 
at angular velocity (2,. with respect to a fixed source, the antenna pattern 
modulation is quite large and polarization dependent via the functions F',(t) 
and F(t) which depend on the orientation of the detector and the source. 

A commonly used parametrization of these amplitude response modula- 
tions is defined in (Jaranowski et al., 1998). 


F,(t) = sin(6) [a(¢) cos(2) + d(#) sin(2¥)] (48) 
F(t) = sin(<) [b(#) cos(2w) — a(t) sin(2¥)], (49) 


where ¢ is the angle between the arms of the interferometer (nearly or precisely 
90 degrees for all major ground-based interferometers), and where w defines 
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the polarization angle of the source wave frame (e.g., angle between neutron 
star spin axis projected onto the plane of the sky and local Cartesian coordi- 
nates aligned with its right ascension and declination directions). The antenna 
pattern functions a(t) and b(t) depend on the position and orientation of the 
interferometer on the Earth’s surface, the source location and sidereal time: 


ae = in(24)(3 —cos(2A))\( — e08(26)) cos2(@— @ — 10,4)1 
-4 cos(27) sin(A)(3 — cos(2d)) sin[2(a — ¢, — 92,t)] 


+5 sin(27) sin(2A) sin(20) cos[(a — ¢, — 2-t] 
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3 sin(27)(3 — cos(2A)) sin(6) sin[2(a — 6, — 92,t)] 
+ cos(2y) cos(A) cos(d) cos[a — dp — 2,t] 


+5 sin(27) sin(2A) cos(6) sin[a — ¢, — 92,t). (51) 


Specifically, in these equations, A is the interferometer’s latitude, and y is 
the counterclockwise angle between the bisector of its arms and the eastward 
direction. The source direction is specified by right ascension a and declination 
6, while ¢,. is a deterministic phase defined implicitly by the interferometer’s 
longitude. These functions reveal amplitude modulations with periods of 1/2 
and 1 sidereal day, and in the case of a(t), a constant term independent of time. 
As a result, the interferometer’s response to a monochromatic source in the 
Earth center’s reference frame will, in general, display five distinct frequency 
components, corresponding to the “carrier” frequency and two pairs of positive 
and negative sidebands, with a splitting between adjacent frequencies of ge ~) 
1.16 x 107° Hz. 

Searches for CW signals must take into account the phase/frequency mod- 
ulations embodied in equation 47 due to detector translational motion and 
the antenna pattern modulations embodied in equations 48-51 due to detector 
orientation changes. 

Figure 9 shows a sample spectrogram of a pure signal simulation using one 
of the so-called “hardware injections” used in LIGO data runs. These signal 
simulations are used to verify end-to-end the detector’s response to a CW sig- 
nal, including sustained phase coherence over long durations. The simulations 
are injected via “photon calibrators,” which are auxiliary lasers shining on mir- 
rors with a modulated intensity. The imposed relative motion of the mirror 
mimics (in the long-wavelength regime) the response of the interferometer to a 
gravitational wave. Various such signals, ranging in nominal frequency from 12 
Hz to 2991 Hz were injected into the LIGO detectors over the O1, O2 and O3 
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observation runs (Biwer et al., 2017). In the example shown, a signal (“Pulsar 
2”) with source frequency 575.163573 Hz (reference time = November 1, 2003 
00:00 UTC) and spin-down —1.37 x 10-18 Hz/s is simulated at a sky location 
of right ascension a = 3.75692884 radians (14h 21m 1.48s) and declination 
5 = 0.060108958 radians (3° 26m 38.36s). Its orientation is defined by incli- 
nation angle . = 2.76 radians (158. deg), and polarization angle = = -0.222 
radians (—12.7deg). The simulation shown (with negligible noise for clarity) in 
Figure 9 applies over a duration of the calendar year 2019 UTC. One can see 
the annual modulation from the Earth’s orbit imposed on an imperceptible 
decrease in the intrinsic frequency. A zoom of 100-hour duration is also shown 
in Figure 10, to indicate the much smaller frequency modulation (by ~2 orders 
of magnitude), along with the intensity modulation. 

As seen from the spectrogram, the frequency modulations lead to station- 
ary bands at the turning points of the modulation. As a result, the spectrum 
averaged over the signal duration peaks at the turning points, as shown in 
Figure 11. These “horns” are a characteristic spectral signature of expected 
signals, where the relative heights of the horns depend on the duration of the 
observation and on the Earth’s orbital phase at the start. For a signal with 
negligible spin-down and a duration equal to a multiple of a year, the horns are 
approximately symmetric, but in the general case that includes observations of 
a few months or less, one or both horns may not be apparent. A more detailed 
analysis of the Fourier transform of a CW signal can be found in (Valluri et al., 
2021). 

In the following, we discuss in more detail the implementations of a selec- 
tion of these methods developed in searches of the initial LIGO and Virgo data 
sets (2001-2011) and that have been further refined for searches of advanced 
detector data. Section 4 presents the results of each type, from searches in the 
data of Advanced LIGO’s and Virgo’s observing runs O1, O2 and O3. 

The volume of parameter space over which to search increases in large 
steps as one progresses through these categories. In each category a star can 
be isolated or binary. Any unknown binary orbital parameters further increase 
the search volume. In all cases we expect (and have now verified from unsuc- 
cessful searches to date) that source strengths are very small. Hence one must 
integrate data over long observation times to have any chance of signal de- 
tection. How much one knows about the source governs the nature of that 
integration. In general, the greater that knowledge, the more computationally 
feasible it is to integrate data coherently (preserving phase information) over 
long observation times, for reasons explained below. 


3.3 Broad approaches in CW searches 


Computational cost depends critically upon the search method used, which 
in turn, depends on the a priori knowledge one has about the source. In the 
following, a broad overview is given of a few key search methods used in 
published searches to date. More details of implementation are presented in 
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Fig. 9 Sample signal spectrogram for a LIGO “hardware injection” (negligible noise), where 
the pixel dimensions are 0.5 hours by 0.556 mHz. 


sections 4, where a selection of these methods is applied to particular classes 
of potential CW sources. 

In this overview, a simplified “toy model” will be used to illustrate scaling 
relations. Methods specific to correcting for modulations will be addressed 
further below in the presentation of particular search implementations. For 
now, amplitude modulation of the signal strength due to rotation of the GW 
detectors with respect to the source is ignored in the following, along with 
Doppler modulations. 


8.3.1 Fully coherent methods 


When applicable, as discussed in section 3.1, fully coherent methods provide 
the best sensitivity. Explicit search methods, taking into account source fre- 
quency and amplitude evolution along with detector noise non-stationarity, 
are discussed in section 4.1. Here, though, let’s consider the highly simplified 
problem of detecting a sinuosidal signal of amplitude hp and known frequency 
fsig, but with unknown phase constant, in random Gaussian noise. 
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Fig. 10 Zoomed-in sample signal spectrogram for a LIGO “hardware injection” (100 hours 
from spectrogram in Figure 9). The sidereal Doppler modulations (~24 hrs) of frequency 
and amplitude modulations (~12 and 24 hrs) are more apparent. 


Imagine the data observation is continuous of duration Tops and has a one- 
sided power spectral noise density function S;,(fa@w), and for convenience, 
assume fig is an integer multiple of 1/Tops (see (Allen et al., 2002) for a 
didactic treatment of the more general case). The DFT bin D; is defined by 
the following sum over a real time series of length Nsample = fsamplefobs With 
values d; (j = 0...Nsample) Where fsample is the sampling frequency of the data 


stream: 
Neample—1 


D; = p> dje71273#/Neample : (52) 
j=0 


From the DFT, one can define the one-sided power spectral noise density 
estimate S),,: 


g,, = MRO + DY") Tae _ 2PM gy 


sample 


for 0 <i < Nsample/2 and where “( )” indicates an expectation value in the 
absence of signal (e.g., determined from an average over many nearby bins and 
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Fig. 11 Sample signal amplitude spectral density for a LIGO “hardware injection” (same 
injection as for spectrograms in Figure 9). 


excluding the bin i itself). Then one can construct a dimensionless detection 
statistic p? using the measured strain power in the DFT bin D; corresponding 
to a signal frequency fsig: 


(54) 


which follows a non-central y? distribution with two degrees of freedom and a 
2 
non-centrality parameter A(ho) = Ayfohe which implies an expectation value 
2 2 7 
2+ crane and variance 4+ 4 Por obs In Gaussian noise one expects a x? distri- 


bution with two degrees of freedom from summing the squares of the normally 
distributed real and imaginary DFT coefficients. 

In the absence of a signal, one can define a threshold value p; corresponding 
to a false alarm probability a such that the cumulative density probability 
function satisfies: 


* 2 


Py 
CDF noise [e7?] = | Parise P; 2) d(p?) =l-a ? (55) 
0 
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where the probability density function is (vy? with two degrees of freedom) 


1 -—2£ 
Dnoise(X; 2) = 5° [2 (56) 


From this threshold and a desired false dismissal rate 6, one can then 
determine the corresponding signal amplitude hae from 


Pi 
CDF signal+noise [px7] = Dupre 2, dh”) d(p?) = B ’ (57) 
0 


where the probability density function is (non-central x? with two degrees of 
freedom) 


1 — 
Dsignal+noise(Z} 2, d) = eae (| Ax) ’ (58) 
and where Jo(y) is a modified Bessel function of the first kind: 
y*/4)? 
Ley SS. (59) 
a 


Choosing a 1% false alarm probability (a = 0.01) leads to p*? ~ 9.21, from 
which numerical evaluation of Equation 57 for a false dismissal probability B 
= 5% leads to an expected sensitivity h2°” of 


Sh, 


Ao” ms 4.544 / 
" Teoh , 


(60) 
which can be taken as a proxy for the expected 95% confidence level upper 
limit on signal amplitude based simply on an observation that an observed p; 
does not exceed p;. In practice, many CW search upper limits are based on 
the loudest statistic found in a search, (e.g., largest p? value for multiple com- 
putations at different frequencies), regardless of whether or not a pre-defined 
threshold has been exceeded. (‘Tenorio et al., 2022) discusses the statistics of 
loudest candidates using extreme value theory. 

The sensitivity expression in Eqn. 60 is shown as the leftmost point of the 
lower blue curve in Figure 12, where it can be compared to sensitivities from 
other methods, discussed below. 

If simultaneous data sets from two independent detectors of identical sen- 
sitivity are added coherently (and a phase correction applied, to account for 
detector separation and source direction), one can construct a combined av- 
eraged detection statistic 7... using the measured power from the square 
of the sum of the simultaneous DFT bin coefficients Di; and Do; containing 


Faiz: 


|Di, ae Dag? Tap 
Di oii = D) js : a — (61) 
NS aes hi 
which has an expectation value of 2+ 2 hoTovs and a variance of 4+ S207 hg Tobe , that 


be 


is, it follows a non-central y? distribution with two degrees of ea and a 
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2 
non-centrality parameter A(ho) = 2g Tee Applying the same methodology as 


above for a single detector, one obtains an expected sensitivity of 


95% — iL S hi 
ho? & Ja x 4.54 Doon’ (62) 
indicating an improvement by \/2 with respect to a single detector. This sen- 
sitivity is shown as the leftmost point of the green curve in Figure 12. More 
generally, N identical detectors with simultaneous data sets’, for which phase 
corrections are known, gain a sensitivity improvement of VN with respect to 
a single detector, as one would naively expect. Put another way, combining 
the N sets gives a sensitivity equal to that of a single detector with an am- 
plitude spectral noise density of \/S;,,/N. Again, this is a simplified model. 
In practice, detectors have different, frequency-dependent 5S), noise levels and 
unequal live times of observing, in addition to different orientations affecting 
antenna pattern sensitivity, ignored here. 


3.3.2 Semi-coherent methods 


Fully coherent methods are computationally costly when covering a large pa- 
rameter space because of the fine steps needed to avoid missing a signal as 
coherence times increase. A crude solution is simply to reduce the coherence 
time and suffer the reduction in strain sensitivity, by approximately Tion. A 
better solution is to apply a semi-coherent method in which the observation 
time is divided into Nseg segments of equal length Teoh = Tops/Nseg, where 
the detection statistic is constructed from an incoherent sum of powers from 
the individual segments. This method sacrifices the constraint of phase con- 
sistency among the different segments and hence is less sensitive than a fully 
coherent method, but is more sensitive than analysis of a single segment alone. 


As shown below, the strain sensitivity scales with 1/Nzig for fixed coherence 
time and large Negcg. 

For illustration, consider a detection statistic constructed from the sum 
of Ngcg individual DFT powers covering the observation time Tops, and once 
again, each normalized by its power spectral density (taken to be stationary 
here, for convenience): 


Neox ) (k))2 
|D; | Teoh 
Rye 4 63 
=1 Ne uh ( ) 


where Neeg = 1 yields R; = p? in Eqn. 54. The underlying statistical distri- 
bution of R; is that of a non-central y? with 2Nsce degrees of freedom and a 


2 
non-centrality parameter of \(ho) = Neeg “42. 


9 Simultaneity is not strictly required, if phase corrections for time offsets between detec- 
tors can be computed precisely enough. 
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In the absence of signal a false alarm probability a implies a threshold RF, 
found from requiring that the cumulative probability density function satisfy: 


R* 


a 


CDF noise [Re] = | Dnoise( Ri; 2) dR; =l-a ) (64) 
0 


where the probability density function is 


gNscg—le—a/2 


= QNecs'(Neeg)’ (65) 


Dnoise(X} 2 Nseg) 
which reduces to Equation 56 for Nscg = 1. RF can be determined numerically 
for arbitrary Nseg, but in the limit of large Nseg, Pnoise reduces to a normal 
distribution with a mean of 2Ngcg and variance 4Nseg, in which approximation 
Ri (a = 0.01) © 2Negeg + 4.654/ Necg- 

In the presence of a signal, the ice value can be obtained numerically 
from the cumulative probability density function: 


RP 
CDF signal+noise [Re] = | Dsignal+noise (Ri; 2, ho) dR; = B, (66) 
0 


where the probability density function is 


~3 (0+ "ye ) S Nseg=* 
signal+ noise x; 2Nse JA = -e 2 Sh ( LOon ) 67 
oa : °) 2 Ned sas ( ) 
LT obs 
INeg—1 (i ou (68) 


which reduces to Equation 58 for Neg = 1, where Tops = NeegZ con has been 
used. 

In the limit of large Nseg and weak signal, however, the distribution ap- 
proaches that of a Gaussian with variance 4Ngcg, from which an approximate 
expression for ae can be obtained: 


Sh 
Tobs ; 


ha ~ VJ/2 [v2 (erfe~ (2a) + erfe“1(28))| ae Nae (69) 


where erfc is the inverse complementary error function. This scaling of sensitiv- 


ity with Nae for fixed observation time is a universal result in semi-coherent 
searches with large Nseg (Krishnan et al., 2004; Mendell and Landry, 2005; 
Prix and Shaltev, 2012; Wette, 2012). It can be understood qualitatively from 
the SNR of an approximately Gaussian detection statistic (large Nscg) scaling 
with ./Neoglcoh = VTovs/Nseg and the direct dependence of that detection 
statistic on the squared signal amplitude. , 

For a = 0.01 and 6 = 0.05, Eqn. 69 yields h95% ~ 2.82N4¢./Sh/Tobs- 
This expression does not agree with Eqn. 60 when Ngcg = 1 because the 
Gaussian approximation breaks down for small Ngcg. For the much lower false 
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Fig. 12 Strain sensitivities (defined by a = 0.01 or a = 10~*; and B = 0.05) of various 
search methods to a bin-centered sinusoidal signal in noise of power spectral density Sp, vs. 
the number of DFT segments into which a fixed observing time Top, is divided. Curves are 
shown for ioe for a = 0.01 and a = 10~7 for semi-coherent searches in a single detector’s 
data (blue) and for a = 0.01 in two-detector searches (green). Asymptotic expressions 
based on the large-Nseg Gaussian approximation are shown as dashed lines. The points at 
the bottom left of each curve represent the fully coherent search sensitivities for 1 and 2 
identical detectors. Semi-coherent curves for two detectors assumed coherent summing of 
simultaneous DFTs for the 2 detectors. Sensitivities for cross-correlation of simultaneous 
data from two detectors are also shown (red). 


alarm probability a = 10-!°, Eqn. 69 yields h95% ~ 4.00N&gV/Sn/Tovs- These 
asymptotic approximations are shown as dashed lines in Figure 12, together 
with numerically evaluated values from Equation 66. 

As is the case for combining data in a fully coherent search from two iden- 
tical detectors with simultaneous data sets, there is a gain of /2 in sensitivity 
for combining simultaneous DFTs from two detectors in a semi-coherent search 
—as long as the DFT coefficients are combined coherently for each segment 
(otherwise, the gain is only 2'/* from semi-coherent combination of DFT pow- 
ers). Figure 12 shows exact (solid green curve) and asymptotic (dashed green 
line) results for two detectors. 

An important variation on semi-coherent methods, used in Hough trans- 
form methods described below, applies thresholding to DFT powers prior to 
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summing. By applying a threshold corresponding to a relatively high false 
alarm rate (and relatively high false dismissal rate for weak signals), one can 
reduce the data volume in processing, to reduce computational cost. In addi- 
tion, by adding integer counts (or, optimally, pre-computed weights) instead 
of measured DFT powers, one’s detection statistic is less susceptible to distor- 
tions from transient, non-Gaussian outliers from instrumental contamination. 
The optimum false alarm probability for thresholding in this idealized sine- 
wave detection analysis (weak-signal limit) is 20% (Allen et al., 2002). 


3.3.8 Cross-correlation methods 


Another attractive approach uses cross correlation between independent (and 
ideally, simultaneous) data streams. The canonical example is cross correla- 
tion between coincident data sets taken with the nearly aligned Hanford and 
Livingston interferometers, but cross correlation can also be used with poorly 
aligned detector pairs and with non-coincident data streams — if sufficient sig- 
nal coherence can be established over longer time scales. In this section, only 
truly coincident data sets will be considered, for simplicity. 

Once again, let’s use the artificial but informative toy model of a bin- 
centered, constant-amplitude sinusoidal signal in Gaussian noise. Let’s also as- 
sume two identically oriented detectors (approximation to Hanford-Livingston, 
which have normal vectors to their planes only 27.3 degrees apart (Althouse 
et al., 1998), in addition to a 90-degree relative rotation about the normal, 
leading to a sign flip in GW response). Also assume that the relative positions 
of the detectors can be accounted for via a signal phase correction for any 
given source direction. In the following, that phase correction is assumed to 
have been applied. 

Using the notation from section 3.3.2, we can combine the two independent 
data streams via DFT coefficients in a narrow band of interest to define a new 
detection statistic: 


2 VaR {D1 D3 ;} Tobs 
N2 Sh, 


sample 


pcc = (70) 


In the absence of a signal, this statistic has an expectation value of zero and 
(by construction) a variance of one. The underlying statistical distribution is 
far from Gaussian, however. In the absence of a signal, one has the sum of 
two normal product distributions (from {D1 iD3,;} =F {D1} R { Do,:} + 
Ss {D1} & {Doi}), which can be obtained analytically, using characteristic 
functions. 

Specifically, a single normal production distribution for the product of two 
zero-mean normal distributions of variance o7 and o3 is 


1 x 
P1 normal product (2) carl Ko ( | | ) ; (71) 


70109 0109 
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where Ko is a modified Bessel function of the second kind, and for which the 
characteristic function is (McNolty, 1973) 


CFL ( le| )i- eee (72) 


Inverting the product of characteristic functions gives the following probability 
distribution for the sum of two such (signal-free) normal product distributions: 


1 |x| 


e€ 7192, (73) 


P2 normal predacs(e) = 
201 02 


For identical detectors (7, = 02 = a) of power spectral density S;,, with 


2 
a signal present of amplitude ho, the expectation value of pac is ee and 


the variance is 1 + ho ng Tobe 

In the presence of a common signal hg in both data streams, one can 
evaluate numerically the value ie for which the false dismissal rate is ( 
for a threshold on the detection statistic corresponding to the signal-free false 
alarm probability a. In the absence of a signal, the threshold on pcc for a false 
alarm probability a = 0.01 is approximately 2.766. For @ = 0.05 one then finds 
Ae = 3.335,/S),,/Tovs, which is only slightly higher than that obtained for a 
fully coherent search using data from two identical detectors (see Figure 12). 

As with semi-coherent searches, discussed in section 3.3.2, one typically 
finds it necessary in wide-parameter searches to segment the data. As before, 
consider dividing the observation time Tops into Nseg equal-duration segments 
of coherence time T.on. In the presence of a signal of amplitude ho, the following 
detection statistic, 


Nseg 


be — 5 Poo: (74) 
Nseg i=1 
hgTcoh _ 1 hoTovs . 1 he aro 
ae. =f vas. and variance Nex [1 + nto 


In the regime of large Neeg and weak signal, the ae probability 
distribution approaches that of a Gaussian for which one expects: 


has a mean value of 


Me) ze, (as) 


nb? w V2 [v2 (erfe-"(2a) + erte-'(28)) "| 2 Tops 


This asymptotic expectation is shown as a dotted red line in Figure 12, to- 
gether with results from numerical simulation over a range of Nseg values (solid 
red curve). This detection statistic is 2'/4 more sensitive than the asymptotic 
1-detector semicoherent sensitivity (lower dashed blue curve), equally sensi- 
tive to the asymptotic 2-detector semicoherent sensitivity for which powers 
from separate detectors are added (not shown), and 2!/4 less sensitive than 
the asymptotic 2-detector semicoherent behavior with coherent summing of 
simultaneous DFTs from the 2 detectors before computing power, as in Equa- 
tion 61 (green dash-dotted curve). 
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One practical consideration to keep in mind for these comparisons is that 
while coherent summing or cross-correlation of simultaneous DFTs provides 
improved sensitivity where possible, those gains are limited by achievable live- 
times of interferometers that operate near their technological limits. 

One can compute nominal signal-to-noise ratios for given signal strengths 
for the coherent, semi-coherent and cross-correlations methods from the noise- 
only variances and the expectation value dependences on signal ho presented 
above. Those SNRs allow sensible direct comparisons among semi-coherent 
and cross-correlation methods when Ngcg is large enough for the noise-only 
detection statistics to exhibit approximate Gaussian behavior over the range 
of interest, but for small Ngcg, including especially the fully coherent case of 
Nseg = 1, the underlying statistics are highly non-Gaussian, requiring care in 
making comparisons. 


3.3.4 Long-lag cross-correlation and loose coherence 


Fully coherent, long integrations seem distinctly different from the multiple- 
short-segment searches based on semi-coherent and cross-correlation described 
(in simplified form) above, and in fact, using fully coherent methods to fol- 
low up on outliers produced by the latter methods is challenging because of 
the typical mismatch in parameter space fineness. Nonetheless, between these 
extremes exist bridges that offer the possibility of both systematic follow-up 
of outliers and more sensitive initial stages for multi-segment searches. The 
first method is long-lag correlation (Dhurandhar ect al., 2008), and the second 
method is known as loose coherence (Dergachev, 2010b). 

Each method benefits from coherently summing DFT coefficients from data 
segments offset in time, which can be motivated by considering the segmenta- 
tion of a single (continuous) Fourier transform of a strain signal h(t) into Neg 
segments: 


F(w; [ta,tz]) = ri a: (78) 
; ee ev iwti-n [iw h(t)e et dt (77) 
Nseg — Teg ti-4a 
Neseg 
1 siwti_ : 
— Nseg 2, € LF (w; [t—15 44) oe) 
1 Neseg 
_ —idi DP. 
=a = e Fi Fi(w), (79) 


where Treg = (tp —t4)/Neeg, ti = 4 +17 ecg and 6; = wt;_1. Hence the Fourier 
transform for the full data span is proportional to the sum of transforms for 
the individual segment transforms with phase corrections e~’®:. 

Now consider once again the artificial but informative special case of a 
monochromatic signal detected via its strength in DFT bins from two detectors 
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1 and 2 with identical observation periods segmented into Nprr epochs for 
which DFT are computed. For simplicity, assume the signal is bin-centered for 
each epoch’s DFT and that the detector noise is both stationary and identical 
for the two detectors. Guided by the relation above, one can then define a 
detection statistic based on the coherent sum of all DFTs from both detectors: 
2 


Norv _ ; 7 ; 
P=| 0 Die + Dae (80) 
t=1 
2 Norr 7 7 : 
= S- S- Dy Di je Hons 999), (81) 
I,J=1 i,j=1 


where $1, and $2, account for the signal phase evolution for each detector 
for each DFT 7 and for any geometric offset between the detectors relative to 
the source direction (see equation 47). This full double sum is computationally 
costly to evaluate explicitly, not only because of the additional operations, but 
more important, because the implicit full coherence requires a fine stepping in 
parameter space. The form, however, makes more clear the relations between 
full coherence and both semi-coherence and cross-correlations (Dhurandhar 
et al., 2008), which can be viewed as subsets of the double sum. A semi- 
coherent sum of powers from individual detectors can be represented by 


2 Noprr 


Psemi—coherent = > s Dr iD}, (82) 


I=1 i=1 


while cross-correlation of simultaneous amplitudes (see equation 74) from the 
two detectors is proportional to 


2 Nprr . _ ; 
Poroae eorélation = S- > Dj 5 ge Moni Oa9), (83) 
UAD)=1 #=1 


This last relation suggests the possibility of following up an interesting 
search outlier from the first stage of a simultaneous-segment cross-correlation 
search by increasing the number of terms kept from the full double sum, al- 
lowing non-simultaneous cross terms and allowing self-correlation terms. For 
example, allowing an offset of up to Njag segment durations would yield: 


2 Nort 
Peross—correlation(Niag) = ye S- Dy gD ge Pr P49), (84) 
IJ=1 45 =1; 
| ~~ j| < Mag 


This approach can also be used at the first stage if rapid frequency evolution of 
the source, such as in a short-period binary system or in young object, argues 
for short DFT coherence times. 

Another, related approach is to define a “loosely coherent” detection statis- 
tic as a subset of the original sum for which phase correlation between nearby 
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(small-lag) DFT coefficients is favored (as opposed to the completely random 
relation allowed by semi-coherent sums). To illustrate’, consider for simplic- 
ity a sum restricted to a single detector. Assume the phase correction applied 
to a product of DFT coefficients for segments separated by a single segment 
lag is taken to be unknown but uniformly distributed in probability between 
—dé and +0. A useful detection statistic can then be formed by (Dergachev, 
2010b) 


1 +6 +6 +6 
Pisese—éoherence _ aren | dA¢i,2 dA¢2,3... BAO Neeg—2,Necg—1 
5 = = 
Necg 
> D; ,D* fe OR aT elgg pate Ae ag!) (85) 
tJ 


where i’(j’) = min(max){i, 7} and A¢gmn = on — dm. Evaluation of the inte- 
grals leads to 


Nseg 


: | 
~ ~, fsin(d 
Piéose coherence = S- Dr, iDz 5 ( \ ’) ’ (86) 
aj 


where the factor (2) ve weights adjacent-lag products more heavily than 
those with longer lags (and yields unity for 1 = j). The single-detector semico- 
herent power sum is recovered by setting 6 = z. In practice, to reduce compu- 
tational cost, the factor is replaced by a discrete kernel function that truncates 
terms for which the weight contribution is too small to warrant the additional 
operations. A generalized version of this detection statistic, including more 
than one detector and non-integer segment lags, has been used in multi-stage 
searches with decreasing 6 at each stage, e.g.,d = 7 > 17/2 > 1/4 > 7/8. 
Each reduction in 6 brings an increase in computational cost per parameter 
space volume as more terms are retained in the sum and the parameter space 
is searched more finely (“zooming in”). 


3.4 Barycentering and coherent signal demodulation 


Taking into account the phase/frequency modulations of equation 47 due to 
detector translational motion and the antenna pattern modulations embod- 
ied in equations 48-51 due to detector orientation changes requires an accu- 
rate model of relevant solar system motions. As noted above, the gravitational 
wave community has adopted techniques of pulsar astronomy researchers, with 
many LIGO data searches using the TEMPO 2 program (Hobbs et al., 2006) 
as a guide and for cross-checking (LIGO Scientific Collaboration, 2018; Abbott 


10 For simplicity, the initial implementation of loose coherence is described here. Later 
refinements (Dergachev, 2012, 2018; Dergachev and Papa, 2019) led to substantial perfor- 
mance improvements. 
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et al., 2004). Correction for the Earth’s and a detector’s motions with respect 
to the solar system barycenter is called barycentering. Independently, Virgo an- 
alysts developed another barycentering software package (Astone et al., 2008), 
also checked against TEMPO 2 and the LIGO software. These packages choose 
steps in time fine enough to allow reliable interpolation of detector motion be- 
tween sampled times. 

Several approaches to incorporating the corrections have been developed 
for continuous gravitational wave searches. These include time-domain het- 
erodyning, Fourier-domain decomposition and hybrid techniques, as discussed 
below. More recently, techniques have been developed for more computation- 
ally efficient barycentering for use in targeted searches (Pitkin et al., 2018) 
and all-sky searches (Sauter et al., 2019). 


3.4.1 Heterodyne method 


Since CW signals are inherently quite narrowband with respect to deviations 
from idealized models, a heterodyning procedure using a base frequency near 
the nominal signal frequency, followed by a low-pass filter allows a large re- 
duction in the number of data samples required to capture the signal modu- 
lations. Conceptually, if one has a pure signal h(t) that can be expressed as 
a slowly varying amplitude function A(t) times a sinusoid of base frequency 
foase, namely, h(t) = R{ A(t)? TfosseT(t)+$0) one can apply the following 
heterodyne for a base frequency fpbase: 


Adi) S hg) e Pn SAG ee, (87) 


where 7(t) relates the SSB time to detector time. This approach allows the 
heterodyned function to have a low effective bandwidth. Applying a low-pass 
filter and then downsampling allows a large reduction in data volume while 
preserving signal fidelity. 

In practice, the heterodyne used in targeted CW searches applies not a pure 
sinusoid factor, but rather a slowly modulated sinusoidal phase ¢yodei(t) de- 
pendent on the topocentric (observatory-centric) time t, a model that includes 
the effects of Eqns. 40 and 47 on signal frequency evolution and propagation 
delays: 

A moaei(t) = d(t) e7 Pmoael (t) (88) 
where it is assumed the signal is well approximated by the model: h(t) = 
R ff (t)elPmoan(t) \, and the data stream d(t) contains h(t) and a (much larger- 
amplitude) random noise n(t). The resulting heterodyne product Hmodei(t) 
can then be interrogated for consistency with noise in addition to a signal 
amplitude function subject to antenna pattern modulations. Small residual 
deviations from the model (“timing noise”) measured empirically from elec- 
tromagnetic pulsation observations ’ can also be taken into account straight- 
forwardly. 

This technique is well suited to searches for known pulsars, for which the 
nominal frequency is precisely known from ephemeris measurements. For ex- 
ample, it has been customary in many LIGO targeted searches to heterodyne, 
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low-pass filter and then downsample to 1 data sample per minute, starting 
from a raw data stream of 16384 Hz. This technique assumes the residual in- 
trinsic bandwidth of the signal following the heterodyne is no greater than the 
Nyquist frequency of 8.3 mHz, which is an excellent approximation for effects 
due to Earth / detector motion. This specific implementation is not as well 
suited to wide-parameter searches, for which the bandwidth must be increased 
or many distinct heterodynes be carried out. 

The resulting heterodyned data samples have had frequency / phase mod- 
ulations due to detector motion removed, but they retain antenna pattern due 
to detector rotation about the Earth’s spin axis. Section 4.1 below presents a 
Bayesian analysis method for such samples (Dupuis and Woan, 2005). 


3.4.2 Resampling methods 


An alternative barycentering technique to heterodyning is to “resample” the 
detector data in order to transform it into SSB time (Schutz, 1991; Jaranowski 
et al., 1998; Brady et al., 1998). A mundane but difficult nuisance is that data 
samples uniformly in detector frame time is not uniformly sampled in SSB 
time, making it difficult to apply conventional discrete Fourier transforms to 
the SSB samples. Two distinct methods have been used to date in CW analysis, 
to make the SSB samples uniform in time. The first (Patel et al., 2010; Meadors 
et al., 2018) uses spline interpolation of the non-uniformly sampled data to 
create uniformly sampled data. In practice, this interpolation is carried out on 
heterodyned subbands, much wider than those used in targeted searches, but 
much narrower than the full bandwidth of the original data collected. Another 
method (Astone et al., 2014b; Singhal et al., 2019) is based on selective data 
sample deletions and duplications, where narrow bands of data are temporarily 
upsampled to much higher frequencies, allowing smaller errors when extra 
samples are deleted or duplicated as SSB time appears to run faster or slower 
than detector-frame time (as defined by successive gravitational wavefronts), 
depending on the relative velocity of the detector with respect to the source. 
As for the heterodyne method, the result in both methods is a data stream 
for which detector translational motion has been corrected, but which still 
contains antenna pattern modulations from daily detector rotation. 


3.4.8 Dirichlet kernel method 


An alternative method can be applied in the Fourier domain by breaking 
the observation time into segments of short-enough duration that the signal 
frequency has negliglible evolution during that duration, that is, the frequency 
change during the time Tcg is small relative to intrinsic frequency resolution of 
a discrete Fourier transform over that duration: Tax" In a templated search for 
a particular signal, the frequency for that segment is known, and a Dirichlet 
filter (Dirichlet, 1829) can be applied to the DFT coefficients in a narrow band 
surrounding the nominal frequency (e.g., +4 DFT bins), using the expected 


weights for those bins for the nominal central frequency. 
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Fractional bin offset of signal from bin center 
0.0 0.1 0.2 0.3 0.4 0.5 

Bin +5: 0 0.0004 | 0.0015 | 0.0030 | 0.0043 | 0.0050 
Bin +4: 0 0.0006 | 0.0024 | 0.0048 | 0.0071 | 0.0083 
Bin +3: 0 0.0012 | 0.0045 | 0.0091 | 0.0136 | 0.0162 
Bin +2: 0 0.0027 | 0.0108 | 0.0229 | 0.0358 | 0.0450 
Bin +1: 0 0.0119 | 0.0547 | 0.1353 | 0.2546 | 0.4053 
Bin 0: 1 0.9675 | 0.8751 | 0.7368 | 0.5728 | 0.4053 
Bin —1: 0 0.0080 | 0.0243 | 0.0392 | 0.0468 | 0.0450 
Bin —2: 0 0.0022 | 0.0072 | 0.0125 | 0.0159 | 0.0162 
Bin —3: 0 0.0010 | 0.0034 | 0.0061 | 0.0079 | 0.0083 
Bin —4: 0 0.0006 | 0.0020 | 0.0036 | 0.0047 | 0.0050 

Sum: | 1.000 | 0.9982 | 0.9933 | 0.9874 | 0.9825 | 0.9807 


Table 1 Fractional powers in neighboring DFT bins (rectangularly windowed) for a 
monochromatic signal with a frequency that ranges from bin-centered (bin 0 of the 10 
bins shown) to a positive offset of a half-bin. The last row gives the total fractional power 
in these 10 bins. 


— a 
[-1Fractional bin offset = 0.0 
(“J Fractional bin offset = 0.1 4 
[-}Fractional bin offset = 0.2 
(“| Fractional bin offset = 0.3 ~ 
{—~] Fractional bin offset = 0.4 _ 
[-}Fractional bin offset = 0.5 | 


Relative Bin Number 


a L a 1 


L L 1 1 L t 
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 
Signal Power fraction 


Fig. 13 Visual representation of the fractional power values listed in Table 1. 


For a bin-centered signal and rectangular windowing, the filter would be 
a Kronecker delta, but in general, spectral leakage favors use of a handful of 
neighboring bins, to recover the full signal strength. By coherently combining 
the resulting extracted complex coefficients from the observed segments, one 
can achieve a full, coherent demodulation of the signal’'. Table 1 shows exam- 
ple power fractions in adjacent DFT bins (using rectangular windowing) for 
a monochromatic signal that is bin-centered or offset positively from the bin 
center by bin fractions of 0.1-0.5 in increments of 0.1. One sees that the bulk 
of signal power can be recovered by a modest number of neighboring bins. 
Figure 13 shows a visual representation of the values in Table 1. 


11 An algorithm commonly used is known as LALDemod (Williams and Schutz, 2000; Prix, 
2018). 
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3.4.4 Time-domain parameter extraction 


After frequency demodulation for detector translational motion, one has a 
highly reduced band-limited, data stream (time-domain for heterodyne or re- 
sampling, Fourier-domain for the Dirichlet filter) implicitly containing the am- 
plitude modulations embodied in equations 48-51, which can be considered de- 
terministic functions of time, dependent upon the putative signal parameters. 
In the case of a source of known position and phase evolution but unknown 
orientation, as for many known pulsars, the unknown source parameters can 
be taken to be the strain amplitude ho, the signal phase constant @¢o, the in- 
clination angle v and the polarization angle w. The data stream can then be 
analyzed to extract those parameters. 

For direct time-domain analysis, the principal method used to date for 
LIGO data analysis has been a Bayesian inference (Abbott et al., 2004; Dupuis 
and Woan, 2005). In brief, the heterodyned data samples {B,} can be ex- 
pressed as complex quantities, with signal template expectations: (adopting 
the notation of (Abbott et al., 2004)): 


y(th,a) = TF ste: th) ho(1 + cos?(z))e'?% 


— 5 Fx (tei th) hha cose, (89) 


where a is a vector with components (ho, 4, t,o) and ty is the time stamp of 
the Ath sample. 

With a set of priors on the a parameters one can extract a joint posterior 
probability density function for these parameters: 


p(al{B.}) x p(a) exp > R{Be — y(tesa)} | 


2 
, 2 OR Be} 


aes Ss S{Br — y(th; a) | (90) 


2 
; 2731 B,} 


where p(a) is the prior on a (uniform for cos(z), ~ and ¢o and ho), and TRL By} 
and as, {B,} are the variances on the real and imaginary parts of By. This 
posterior distribution can be examined for evidence of a signal present. In the 
absence of a signal, an upper limit on strain amplitude hg can be found via 
marginalization over the other three signal parameters to obtain a marginalized 
posterior: 


p(hal{Bu}) « II p(al{Bu}) de dy) do, (91) 


normalized so that i, p(ho|{Bx})dho = 1. Unlike a frequentist confidence 
level, the resulting curve vs ho represents the distribution of degree of belief 
in any particular value of ho, given the signal model, the parameter priors and 
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the data observations {B,}. One can derive a 95% credible Bayesian upper 
limit 85% U" for which the probability lies below h9°% UE via 


n25% UL 


0.95 = | p(hol{Be})dho. (92) 


The combined posterior distribution from multiple, independent detectors 
can be obtained via the product of the individual likelihoods (Dupuis and 
Woan, 2005). In the event that estimates of » and 7 can be inferred from 
electromagnetic measurements of the source, e.g., from images of jets assumed 
to be emitted along the spin axis of a star, then the precision on hg can be 
improved by assigning much narrower priors to the parameters. 


3.4.5 Five-vector method 


The so-called “Five-vector” method exploits the property that the complexity 
in equations 46 and 48-51 can be distilled down to five terms (Astone et al., 
2010b) 

h(t) = hoA - We'iot+9o) | (93) 


where wo is the signal frequency in the SSB frame, where A can be decomposed 
into plus- and cross-polarized terms that depend on complex amplitudes H+ 
and Hy: 

A=H,At+H,A%, (94) 


and where At and A* can be expressed in terms of trigonometic functions, 
using equations 48-51 (see (Astone et al., 2010b) for detailed expressions). 
The vector W is a five-component set of basis functions, indexed by k = 
{—2, —1,0, 1, 2]: 

W, =e 9, (95) 


where OQ is the detector’s local sidereal time in radians. 
The data stream x(t) too can be decomposed using these basis functions: 


x= | a(t wert (96) 


One can then construct a detection statistic using a weighted sum of the 
squared projections: 


S=cy|hy? +ex|ax/?, (97) 
where the projections are defined by 
: X- At 5 X- At 
pi At?’ or Ax?” (98) 


Empirically (Astone et al., 2010b), it is found that best performance for known 
1, w can be obtained with the weightings: c,,. = |At™|*, while estimation of 
signal amplitude can be obtained from 


ho = V/|h4l? + |x. (99) 
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3.4.6 The F-statistic 


The most pervasive detection statistic used in broadband CW searches can 
also be used for targeted searches, namely the F-statistic (Jaranowski et al., 
1998). As above, the F-statistic is constructed to take into account not only 
the frequency / phase modulation of the detector’s translational motion (using 
time-domain or frequency-domain techiques), but also the amplitude modula- 
tion from daily detector rotation. 

It is constructed from a general maximum likelihood approach, where the 
data is taken to be a sum of random noise n(t) and a signal h(t): 


a(t) = n(t) + A(t), (100) 
where h(t) from equations 48-51 can be written!” 

h(t) = >> Ajhi(t), (101) 
where the coefficients A; are inferred from equations 48-51: 


A, = hosin(¢) stl + cos*(z) cos(2~) cos(2 By) 


—cos(c) sin(2w) sn(2%) : (102) 


Ag = ho sin(¢) (l + cos*(z) sin(2w) cos(2 p) 


+ cos(z) cos(2w) sn(20) : (103) 


Ag = hosin(¢) = + cos?(z) cos(2w) sin(2 Bo) 


—cos(c) sin(2w) cn) ; (104) 


Ag = hosin(¢) F(t + cos?(v) sin(2w) sin(2 By) 


+ cos(z) cos(2w) cn) ; (105) 


12 In the original F-statistic article (Jaranowski et al., 1998), a two-component signal model 
is assumed, corresponding to frequencies at once and twice the source rotation frequency. 
Only the component at twice the rotational frequency is considered here where the wobble 
angle @ in (Jaranowski et al., 1998) is taken to be 7/2 for a triaxial ellipsoid, allowing a 
simplification of notation. 
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and the time-dependent functions h; have the form: 


hy = a(t) cos(2(t)), ——-Aa(t) = b(t) cos(20(t)) (106) 
hg = a(t)sin(20(t)), -hg(t) = b(t) sin(2H(t)), (107) 


where &(t) is the phase of the signal, including modulations. 
A log-likelihood function log(A) is constructed via: 


log(A) = (21h) — 5(Alh), (108) 


where the scalar product (|) is defined by a filter matched to the detection 


noise spectrum: Sor KS 
(2ly) = anf 7 Or a}, ue) 


where ~ denotes a Fourier transform, * is the complex conjugation, and Sj, is 
the one-sided power spectral density. 

Following (Jaranowski et al., 1998), the narrowband signal allows, in prin- 
ciple, conversion of the scalar product to a time-domain expression: 


9 Tobs 
(ais smi x(t)h(t)dt, (110) 
Sn(fo) Jo 

where stationarity of the noise over the observation period Tops is implicitly 
assumed, which unfortunately, is rarely a good assumption for interferometers 
at the frontier of technology. Nonetheless, practical implementations of the 
F-statistic are not limited by this assumption. Defining a time-domain scalar 

product: 


) Tobs 
(ally) = p— [a ute, (111) 
obs JO 
the log-likelihood function can be approximated via 
Tobs | 1 | 
log(A) = x\|h) — =(h\|h)] , 112 
e(A) Sth) | 5 (hi|h) (112) 


which is proportional to a normalized log-likelihood log(A’): 
1 
log(A’) = (a||h) — 5(Allh), (113) 


which does not depend explicitly on the spectral noise density. The signal 
depends linearly on the four amplitudes A; and can, in principle, be extracted 
from a likelihood maximization: 


Olog A’ 
5 7 =0, (114) 
from which a set of linear algebraic equations can be derived: 
4 
> Mi A; = (2A), (115) 


i=1 
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where the components of the matrix M;; are given by 
Mi; = (hillhy). (116) 


Cross-terms of the cos(®(t)) and sin(®(t)) terms in equations 106-107 can be 
neglected in the time integrations. The surviving terms can be expressed: 


(hy|h1) © (hglh3) & =A, 
(ha|h2) © (hilha) © SB, 
(hilh2) & (hg|ha) © =C, (117) 


where A := (alla), B := (b||b) and C := (al|b). With these approximations, 
the matrix M becomes 


CO 
Meloe (118) 
where O is a zero 2 x 2 matrix, and C is 
1/AC 
c=5(Ge) (119) 


from which maximum-likelihood estimators A; of the true amplitudes A; can 
be obtained: 


Ble hy) = C(a hg) 


A = 7 ; 
Ay = 9A ba) = Cel) 
Ag = 2B ta) = Calf) 
Ay = 2A Vie hs) et) 


where D = AB — C?. Substituting these expressions into equation 112 leads 
to the F-statistic (denoted by 2F: 


B(a||h1)? + A(a||h2)? — 2C(a||h1)(a|Ih2) 
= Bn (fo) D 


, Blellhs)? + Aaa)” = 2C(2||hs)(2/|ha) 


_ (121) 


The quantity 27 has a probability distribution of a chi-squared with four 
degrees of freedom in the absence of a signal and that of a non-central chi- 
squared with a non-centrality parameter: 


d= d* = (Alh) (122) 
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where d is proportional to signal amplitude (Jaranowski et al., 1998). The 
probability distributions Ppoise(2F) and Psignal+noise(2F; d) are hence: 


1 
Pnoise(2F) = {OF OPM, (123) 
2 2 1 172 
Dermnali nce ee d) = ( a Ty (ay (@F)) eae i-ae ) (124) 
(125) 


where J; is a modified Bessel function of the first kind (order 1). 
As discussed in section 3.3, 27 can be used as a detection statistic, where 
a threshold 27 can be chosen to satisfy a desired false alarm probability: 


2Fo 
CDF noise[2Fo] — | Pnoise(2F) d(2F) =1- Qa, (126) 
0 
1 1 
=1- (1 +2Fo+ 5270 + 728) arn, (127) 


and where the probability of detection for a given d is 


Co 


Paetection(d, 2Fo) = | Dsignal+noise(2F; d) d(2F). (128) 
2Fo 


The formalism above describes a time-domain implementation (Jaranowski 
et al., 1998; Astone et al., 2010a), but a narrowband frequency implementa- 
tion (Prix, 2018) has been used extensively in LIGO searches. 

In searches for known pulsars for which optical or X-ray observations of 
pulsar wind nebulae allow inference of 4 and w, a modified version of the F- 
statistic known as the G-statistic can be applied to gain slightly in sensitivity, 
depending on the stellar orientation (Jaranowski and Krolak, 2010). 

Although originally derived in a frequentist, log-likelihood framework, the 
F-statistic can also be obtained in a Bayesian approach (Prix and Krishnan, 
2009) with an unphysical prior (non-isotropic in stellar orientation), an alter- 
native framework that has received additional study (Prix ct al., 2011; Keitel 
et al., 2014; Whelan et al., 2014; Dhurandhar et al., 2017; Bero and Whelan, 
2019; Wette, 2021). 


3.5 Semi-coherent signal demodulation 


Let’s now consider a coarser demodulation, in which phase fidelity is not re- 
quired for the full observation time. Instead, the observation is broken into 
discrete segments of coherence time T,., which need not be contiguous with 
each other. The segmentation reduces the fineness with which the parameter 
space (e.g., frequency, frequency derivatives, sky location) must be sampled, 
leading to often dramatic reduction in computing cost to search a given param- 
eter space volume, albeit with a degradation of achievable strain sensitivity. 
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df/dt 


Frequency 


Fig. 14 Conceptual illustration of the “stack-slide” method in which rows of a spectrogram 
are shifted up or down in frequency to account for Doppler modulations. 


3.5.1 The stack-slide method 


For short-enough T,on, no frequency demodulation need be applied within a 
single segment. One can simply sum the strain power from each bin in an 
DFT containing the frequency of the signal for that time interval. Figure 14 
illustrates the simplest version of this approach, known as “stack-slide” (Brady 
et al., 1998; Mendell and Landry, 2005). In a spectrogram where each column 
represents DFT powers for a given Tyon, the bin containing the signal fre- 
quency (indicated by the green square) varies in frequency from one column 
to the next. Correcting for the frequency modulation by shifting columns up 
or down leads to the signal’s power being contained in a horizontal track in 
the demodulated spectrogram. For a relatively narrow frequency band, the 
amount of vertical shift for a given column is nearly the same for all frequen- 
cies in the band, for a given set of source parameters, including sky location. 
Hence by stacking powers across rows in the demodulated spectrogram, one 
can look for an outlier indicating a signal. 


To be concrete, define the power p®) to be the strain power spectral 


density measured in bin i of DFT k, where the bin i is the appropriate bin 
after “sliding”: 


ry(*) 
2|Di(demoa) |? 


p(k) 
P. — 
: Teoh 


(129) 
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Following (Abbott et al., 2008a; Mendell and Landry, 2005), this power is 
renormalized to form a dimensionless quantity 7* 


p®) 


= 1 130 
=a (130) 


ni 


where S$ h is the one-sided power spectral density expected in the absence of 
signal. This quantity differs from the p? defined in Equation 54, both in the 
implicit demodulation associated with bin i and in a factor of 2. Here 7* has 
an expectation value of 1 in the absence of signal. 

The stack-slide detection statistic Pi then is the average value of 
over the Nprr DFT’s used in the analysis: 


(e) 1 Nprr 
Prey = ’ 131 
i(SS) Nore », UT ( 3 ) 


This quantity has an expectation value of 1 in the absence of signal and a vari- 
ance of 1/Nprr. Signal candidates are chosen based on exceeding a threshold 
corresponding to a false alarm probability a, from which detection sensitivity 
is determined from a desired false dismissal probabiity 8. Appendix B of (Ab- 
bott et al., 2008a) details the statistical behavior. In brief, the quantity (similar 
to R; of equation 63 above) 


Piss) = 2Nprrnt (132) 


has the probability density distribution of a non-central x? with 2Nprr de- 
grees of freedom and a non-centrality parameter 2Nprr <d?> which is the 
expectation value of the estimator in Equation 122 when evaluated over a 
single DFT. Hence the probability density distribution for P;(gg) follows: 


INppr—1 (\/Piss)Nprr <2 >) 
(Nprr <d? See 


Nprr-t d? 
N, <-> 
z pFTt+<% 
x Piss) € ( ), 


Psignal+noise (Piss); Nort; d) > 


(133) 


Numerical evaluation (Abbott ct al., 2008a) for a = 0.01 and 8 = 0.10 leads (in 
the large Nprr limit) to a sensitivity <d? >@ 7.385/./Nprr and to a strain 
sensitivity for a single template search of pi) se 7.74/S),/ any Be a where 
Tops refers here to the total observing time analyzed and where stationary 
data is implicitly assumed. In practice, however, this method is applied to 
wide-parameter searches for which trials factors lead to much worse strain 
sensitivities (Tenorio et al., 2022). (Prix and Shaltev, 2012) carry out a detailed 
analysis of maximizing sensitivity at fixed computational cost for different 
stack-slide search configurations. 
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8.5.2 The Powerfluz method 


The PowerFlux method (Abbott et al., 2008a), in its simplest form, is similar 
to the stack-slide method, with the following refinements: 1) an explicit polar- 
ization is assumed for each signal template searched, with an antenna pattern 
correction applied; 2) detection statistic variance is minimized in the presence 
of non-stationary noise; and the detection statistic itself is a direct measure of 
strain amplitude. 

Using the same notation as above (see Equation 129), the PowerF lux detec- 
tion statistic Rpp for a given set of orientation parameters z and w is written: 


2 here wiPM” (Fie, ¥))? 


Rpp= 134 
PP= 7p rr yy, ; (134) 

where the weights are defined as 
Wi = (File, d))*17/ Sh, (135) 


and where F;(v,w) is the antenna pattern weight calculated for the midpoint 
of the time segment 7 for the assumed polarization such that the detector 
amplitude response can be written as haet,; = hoFi(e, y). In practice, searches 
have been carried out for circular polarization (t = 0 or 7) and for particular 
linear polarization angles w (with 1 = 7/2) to define “best-case” and “worst- 
case” orientations, respectively. 

The choice of weight definition comes from minimizing the variance of the 
strain amplitude estimator p® /(Fi(u,v))?, where the noise (in the weak sig- 
nal regime) is assumed to be dominated in each time segment i by a power 
spectral density S;,, with underlying Gaussian distributions for real and imag- 
inary DFT components. Under that assumption, the variance of the noise is 


proportional to (S;,,)?. As a result, each term in the numerator of Equation 134 
is proportional to (F;(u, v))?P™ /S;,,, which gives higher weight to segments 
with higher F(z, ~) magnitude and lower noise 5;,, as one would wish. For a 
given polarization choice defined by (z, 7) the detection statistic Rpp is a di- 
rect measure of total strain power such that subtracting the expectation value 


based on neighboring bin yields a direct estimator for signal power. 


3.5.3 Hough transform methods 


Hough transform methods refer, in practice, to an application of a pattern 
recognition algorithm first developed for use in the 1960’s by high energy 
particle physicists (Hough, 1959, 1962) to reconstruct a charged particle’s tra- 
jectory from discrete positions (“hits”), measured by a tracking detector. The 
method is best suited to data that is “sparse” and for which a simple transfor- 
mation from the raw measurements to the signal parameter space can amplify 
the detection statistic. In the original application to particle tracking, the 
hits were two-dimensional projections for which looking for straight lines built 
out of all hit combinations was computationally intensive (especially in the 
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1960’s!). To represent a straight line, instead of offset and slope, the vector of 
its minimum distance to the origin, in polar coordinates (r, #), is used. A point 
(x,y) belonging to that line sets the relation r = xcos@ + ysin@ which is a 
sinusoidal curve in the @-r plane. Cells in that plane count how many curves 
pass within their boundaries, and the most occupied cell identifies (r, @) of the 
original track. 

In the case of CW searches, two different Hough transform methods (“Sky 
Hough” and “Frequency Hough”) have been used in recent years, both of 
which accumulate excess power from frequency-demodulated DFTs. In the Sky 
Hough method (Krishnan et al., 2004), the transformation is from a narrow 
frequency band and frequency derivative to right ascension and declination, 
where broad patches of sky are searched collectively. In the Frequency Hough 
method (Antonucci et al., 2008; Astone et al., 2014a), the transform is from a 
time-frequency plane to a plane of frequency and frequency derivative. In each 
case, one searches for a statistically significant excess among the pixels and 
applies a thresholding to individual accumulated powers, in order to reduce 
computational cost in the accumulation. 

The Hough number count is defined as a weighted sum of binary counts 
4: 


n= S- WiNi, (136) 


where n; = 1 if the normalized segment power 7* exceeds a threshold 7* and 
zero otherwise, and where the weights favor low-noise times and are optimized 
for circular polarization (Antonucci et al., 2008; Abbott et al., 2008a): 


wi x s— (Et)? + (FX)? ], (137) 


with a normalization chosen to satisfy: 


Nprr 


S- Wi= Nort. (138) 
i=1 


In the Sky Hough method (Krishnan et al., 2004), so-called “Hough maps” 
in right ascension and declination are created for each assumed frequency and 
frequency derivative, where signal outliers produce “hot” pixels in the sky 
patch for which the map applies. In the Frequency Hough method (Antonucci 
et al., 2008; Astone et al., 2014a), the Hough map is created instead in the plane 
of frequency and frequency derivative for each localized sky point. The primary 
motivations for this alternative mapping to parameter space are reduction of 
inaccuracies arising from approximations and non-linearities in the mapping 
to the sky and avoidance of artifact “pileup” in which certain regions of the 
sky are contaminated over subbands by particular narrowband artifacts. 

Regardless of the choice of parameter space mapping, the statistical char- 
acter of the Hough number counts is governed by the value of the threshold 
used to define the binary counts n;. The mean number count in the absence of 
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a signal is mn = Nprrp, where p is the probability that the normalized power 
nk exceeds a threshold value n*. For unity weighting, the standard deviation 
is on = Nprrp(1 — p). For the more general weighting, this becomes: 


Nprr i 
| (139) 


o= pu» Yo wi 


i=l 


For Nprr > 1, the underlying distribution can be approximated as Gaussian, 
in which case a threshold n,(a@) corresponding to a false alarm rate a is given 
by (Krishnan et al., 2004) 


Nth = Nprr p+ V2oerfc' (2a), (140) 


where it is natural to regard the significance of a given measured n to be 


(141) 


In (Krishnan et al., 2004; Abbott et al., 2005a, 2008a) an optimal choice of 
the normalized power threshold parameter is found to be n* = 1.6, for which 
p=ew” = 0.2. 

One can compute (Abbott et al., 2008a) a sensitivity ha (a) for a false 
dismissal probability 6 and false alarm probability a: 


n}-8(a) wv 3.38(8)1/2 ( hwll a a (142) 
° , w-X Tots 
where ||w|| = S72"? w? and 
S = erfe~!(2a) + erfe™ (28) (143) 
1 4 2 a 2 
X= 5 (Fi)? + (FY), (144) 


and where He Ise refer to the antenna pattern functions for the + and x po- 
larizations evaluated at the midpoint of time segment 7. 

Other improvements to the Sky Hough method have included incorporating 
a hierarchical approach (Sancho de la Jordana, 2010), adaptation to a search 
for stars in binary systems (Covas and Sintes, 2019) (see section 4.5), clustering 
of outliers (Tenorio et al., 2021c) and systematic outlier follow-up (Tenorio 
et al., 2021a). 


8.5.4 The Stacked F-statistic method 


The semi-coherent approach used above (in various approaches) with DFT co- 
efficients can also be applied to longer segments of time for each of which the 
coherent F-statistic is computed. This approach permits deeper sensitivity 
since the F-statistic can be computed without degradation of signal coher- 
ence for arbitrarily long periods of time. The disadvantage is that the much 


72 Keith Riles 


finer resolution in parameter space associated with such sensitivity leads to 
much greater computational cost, coming from the fine stepping needed within 
each segment and from the mapping with negligible signal loss from one seg- 
ment to the next. A variety of F-statistic “stacking” methods'’ have been 
implemented over the years, both inside and outside of the framework of the 
Einstein@Home distributed computing system (see section 4.4). When com- 
puting the F-statistic over short time segments, a modified variation, the Fap- 
statistic, which avoids degeneracy due to minimal antenna pattern modulation 
can be more effective (Covas and Prix, 2022). 

Many of the considerations discussed in semi-coherent summing of DFT 
power have analogs in F-statistic summing, including the use of threshold- 
ing and the use of Hough transform mapping. Particular implementations will 
be discussed below in sections 3.6.1 and 4.4. One critical issue in these com- 
putationally costly searches is the optimum placement of signal templates in 
parameter space, to be discussed next, more generally. 


3.6 Template placement 


Computationally demanding searches must choose step sizes in signal parame- 
ter space, with finer spacing leading to greater cost, in general. The choices are 
typically governed by what is considered an acceptable maximum “mismatch”, 
normally parametrized by the fractional decrease in detection statistic for a 
given offset in parameter space. 

For an n-dimensional, locally Cartesian grid defined by n search param- 
eters, one can regard the mismatch parameter js as governing the maximum 
half-length of the diagonal of the n-dimensional cell containing the correct 
signal parameters. Conceptually, we imagine having made the least optimum 
choice of grid offset such that the true parameters lie at the center of the cell, 
and no matter which of the 2” corners of the cell is sampled, the value of the 
detection statistic is no smaller than 1 — yu of the value obtained, had the cen- 
ter of the cell been sampled. Figure 15 illustrates the concept with a detection 
statistic “surface” above a plane in two signal parameters, where the contours 
correspond to mismatch values of 20%, 40%, 60% and 80%. 

In the following, general considerations of template placement are consid- 
ered, first for directed searches for particular points on the sky, for which 
placement is relatively straightforward, and then for all-sky searches, where 
template placement is quite subtle and remains an active research front. 


8.6.1 Template placement in directed searches 


For coherent directed searches, the phase evolution equation 40 governs tem- 
plate placement, where for multi-day analyses, the effects of amplitude mod- 


13 “Stacking” the F-statistic values is more subtle than in the stacking used in the stack- 
slide and other semi-coherent methods based on summing DFT powers because the demod- 
ulations to obtain the F-statistic values differ across time segments. 
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Fig. 15 Illustration of mismatch for a generic detection statistic. The upper panel shows a 
“surface” of height equal to the detection statistic for a pure signal above a plane defined by 
two signal-defining parameters (with zero covariance for simplicity). The green cross marks 
the true location for the two parameters and the maximum possible detection statistic. 
The lower panel shows detection statistic contours in the two-parameter space, where the 
contours correspond to mismatch values of 20%, 40%, 60% and 80%. The red crosses define a 
search template grid chosen to be least optimal for this signal location in that the true signal 
location is centered in a 2-dimensional cell, which maximizes the possible minimum mismatch 
(20%) between the detection statistics for the true signal and the closest template. The 
dashed diagonal line defines the “distance” in the 2-dimensional parameter space between 
the true signal location and the closest search template. 


ulation can be safely neglected in choosing template spacing (Prix, 2007a,b). 
Consider for a moment a highly simplified detection statistic based on mul- 
tiplying in the time domain an assumed sinusoidal signal template having a 
particular phase constant @% and frequency fj against the raw data x(t), as- 
sumed to be a sum of random Gaussian noise n(¢) and a sinusoid signal having 
amplitude ho, phase constant ¢9 and frequency fo: 


2 
PP aay . 
Fast) = 7 | e F042 Fol) a(t) dt (145) 
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T 
nm | e— (042708) [n(t) + ho cos(¢o + 2m fot)|dt| . (146) 


0) 


In the limit of large T’ and strong signal, the expectation value of F when 
maximized over possible template values for ¢j and f{ is simply h2, inde- 
pendent of $9 and fp, where F is maximized for Ad = ¢’ — ¢ = O and 
Af = f$—fo = 0. To understand how rapidly F’ decreases as | A f| departs from 
zero, it’s helpful to rewrite cos(¢o + 2m fot) = 5 (e%(Got?t fot) 4 etldo+2n fot) ) 
where in the strong-signal limit of large T, F approaches 


ie 2 
p= |" | eer dt (147) 
T Jo 
2 
= h2 | _ ge trar) (148) 
9 | (Qa AfT 
2 2sin(2rAfT/2) |? 
_ 72 ad 14 
ho le | In AFT a 
= h2 |sine(n AfT)| (150) 
1 
= h2 E ae (rajr)’| ; (151) 
where the convention sinc(#) = sin(a) is chosen. If we rewrite this last result 


as F = h3cos*(A¢@mismatch), then the tolerance in Af for a phase mismatch 
value Admismatch is 


V3 
Af mismatch Ly Admismatch; (152) 
aT 


which is 2\/3 larger than the naive underestimate of equation 41. Consequently, 
in a search that automatically maximizes F’ over the unknown phase constant, 
one need not search as finely in frequency as suggested by equation 40, which 
implies reduced computational costs in large-scale searches. 

Given the importance of template placement to those costs, in fact, a sys- 
tematic approach is merited. Following methodology developed originally for 
template placement in compact binary merger searches (Sathyaprakash and 
Dhurandhar, 1991; Owen, 1996; Balasubramanian et al., 1996), one can rewrite 
and generalize the simplified detection statistic in equation 145, replacing the 
data with another template and address the reduction in F’s value due to 
mismatch of template parameters 


2 
if T j AY) 5 7 
Plo, 83) =| fo ewiPOA DCO at (153) 
0 
1 2. ' ° 
_[i | ei(P(EA)—G(A')) gy (154) 
T Jo 
iy ae ; 
= zi ei AP(HAA) Gy (155) 
T Jo 
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where A and X’ refer to a set of N parameters, such as phase and frequency 
derivatives, and where AA = A — X’. Clearly, for AA = 0, F = 1 and is 
maximum, with vanishing first partial derivatives. Hence we expect F' to have 
the following form in the vicinity of AXA = 0: 


Fei Ad; AA 156 
> ner Oe i ey 
which leads to the definition of a metric: 
1 OF 
.. = —-——__—_ 157 
ote 2 OAM OA | nse cae 


such that the mismatch pw of a template deviation is u = 3°) » greAAnrAre. 
Hence the appropriate spacing of templates in parameter space to avoid ex- 
cessive mismatch is governed by the form of gxe. 

A general treatment of finding gzg (Owen, 1996) can be approached by 
Taylor-expanding the exponential in equation 153: e'4? ~ 1+ iA’ — 5 AG? 
and evaluating the second derivatives of F' with respect to AA, and Adg. In 
the limit AX > 0, one finds: 


_ | a ae aD af Aids 
Arco \OA% BAN! \8ADg/ \OAX2/? 


TE 
(=z fp Hoe ieee ae 


More specifically, in the context of the Taylor expansion of the phase function: 


OF 
OAX,OAX¢e 


where 


(159) 
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where f(™) = FA |) o F can be expanded: 
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m,n=0 


Terms higher in order than Af(™ Af) have been neglected in the above. 
From this last expression, we conclude that the metric gze can be written: 


Tete+2(k + 1)(€+1) 
(K+ 2) 0+ 2)'(k +643) 
See (Wette et al., 2008) for the same expression for the metric for the F- 
statistic (Jaranowski ct al., 1998) in a directed search. 


As examples, consider the 2-parameter metric with respect to frequency fo 
and its first derivative f;: 


Gke = (27)? (166) 


1 

900 = 3 (xT)? (167) 
1 2 

901 = & (ar/?) (168) 
4 

Hi = 3, (nT?)" (169) 


For a given desired mismatch AM, define nominal offsets Afj and Afy, using 
only the diagonal metric elements: (Aff = VAM /gxx) 
V3AM 


Afi = (170) 


»  38V5AM 

Af = (171) 
Since off-diagonal terms in the metric are non-zero, a rectangular grid using 
only diagonal terms will, in general, be inefficient. Figure 16 illustrates for a 
2-dimensional slice of Afp vs. Afi (=A faw) a template grid that accounts for 
these correlations in mismatch. A grid placement based on only the diagonal 
metric elements would lead to inefficient coverage, as shown. References (Prix, 
2007b; Wette, 2014) discuss more generally and in more detail template grid 
placement for CW searches, with special focus on searches over the three- 
dimensional parameter space (faw, few, faw). As noted above, however, for 
short coherence times, the range of few searches may be smaller than A eon 
in regions of parameter space, depending on braking-index assumptions and 

the value of faw. 
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Fig. 16 Illustration of a (Afo, Afi) template grid (black stars) and constant-mismatch 
elliptical contours for which the grid point placement gives complete coverage. The values 
of the frequency and frequency derivative are given in normalized units of Afj and Aff 
defined in equations 170-171. The magenta diamonds indicate a rectangular grid with full 
coverage when the off-diagonal metric term is ignored. 


8.6.2 Template placement in all-sky searches 


Template placement in all-sky searches is relatively straightforward for semi- 
coherent searches using short coherence times Ton of ~(hours) or less, but is 
quite subtle for coherent searches using much longer coherence times and for 
semi-coherent searches using long coherence times for each data segment. 
Short-T.on template grids can be factorized over sky location (a, 6) and over 
(faw, faw), using isotropic grid point placement, e.g., density proportional to 
cos(6) and uniform in a, with a rectangular grid in (few, few), with spacings 
determined empirically or semi-analytically for a given data run. For example, 
the rule of thumb given in equation 44 overestimates the density needed for 
short observation times of ~(few months) because of correlations (Prix and 
Itoh, 2005) in the dependence of a semi-coherent power sum on sky location 
and frequency parameters. For a data set collected over 1-2 months of the 
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Earth’s orbit, the average acceleration of the detector toward the Sun creates 
an apparent offset in the spin-down of a putative source. Hence a search over 
a band of frequencies and 1st derivatives may detect a signal with nearly as 
high an SNR as the nominal maximum, but with correlated offsets in the 
four parameters (a, 6, faw, faw). For longer observation times, these near 
degeneracies in parameter space become less helpful; signal templates must be 
placed more densely. At additional computational cost, these semi-coherent 
searches may also search explicitly over source polarizations, or may choose 
to apply a circular polarization weighting and sacrifice some sensitivity to 
near-linear polarizations (Abbott et al., 2008a). 

Template placement for much longer coherence times is more challenging 
because analytic approximations break down for long coherence times and 
because naive grid spacings depend on the specific region of the Earth’s orbit 
covered by a particular coherence time, making the systematic matching of 
signal candidates across different time segments non-trivial in semi-coherent 
searches. Template placement for the F-statistic has received much attention 
in the last decade and a half (Whitbeck, 2006; Prix, 2007a,b; Wette and Prix, 
2013; Wette, 2014), in part because of its use in the Einstein@Home (see 
section 4.4) distributed computing platform. From equation 122, one can define 
a mismatch analogous to that of section 3.6.1: 


— (hayaalha+aa) 
—  (AylRy), 


where we expect an approximately quadratic falloff from unity for small | AA| 
(but see (Allen, 2019) for a discussion of template placement for larger |AA| 
and see (Allen, 2021) which distinguishes between optimality for setting rig- 
orous upper limits and optimality for signal detection). 

The complexity of the definition of (h|h) (see equations 50-51 and 101- 
107) do not yield a definition of (h|h) in the convenient form of equation 153. 
In particular, the sidereal antenna pattern modulations due to the Earth’s 
rotation are not accommodated by the phase-only dependence of the simplified 
form. For long observation times, however, amplitude modulation effects can 
be averaged with sufficient accuracy (Prix, 2007a). Phase modulation from the 
Earth’s motion is captured by equation 153, allowing use of equation 158 to 
determine the F-statistic metric with respect to frequency parameters and sky 
location. 

Following the treatment of (Prix, 2007a), a more explicit phase evolution 
can be written: 
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(172) 
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N p(m) 
P(t) = Go +20 > i ; (173) 
m=0 


where r(t) is the SSB arrival time of the signal. Ignoring the Shapiro and 
Einstein delays in equation 47 for metric definition, one can write: 


r(t)- A 


— Tref ; ( 1 74) 
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where r(¢) is the position of the detector at time t, * is the unit vector pointing 
from the detector to the source, and Tye is the reference time in the SSB frame 
at which the frequency and its derivatives are defined. 

The phase derivatives entering equation 158 can then be written: 


ab r(t)h+ 

dF® ~ “"(k+ I) (175) 
Ob rift) Cy fO™ (ter) (t)™* 
ni oe c x (m+ 1)! , ites) 


from which the metric terms for a particular point in parameter space (f, 
f) can be computed via numerical integration of equation 158 over the ob- 
servation span, with precise description of r(t), accounting for the non-zero 
eccentricity of the Earth’s orbit. It is convenient in some studies, though, to 
make the “Ptolemaic” approximation (Jones et al., 2005; Whitbeck, 2006) in 
which the Earth’s orbit is treated as circular, for which analytic but quite 
lengthy trigonometric expressions can be obtained (Whitbeck, 2006). 

An inconvenient property of the metric defined above using the param- 
eters (f, #) is that converting the 3-D Cartesian # components to the 2-D 
sky coordinates a and 6 leads to a sky spacing that depends on the param- 
eter themselves. A metric more convenient for large-scale CW searches over 
the entire sky can be obtained by exploiting global correlations in parameter 
space (Pletsch, 2008; Pletsch and Allen, 2009). For multi-day coherence times 
short compared to one orbital year, one can Taylor-expand the rescaled posi- 
tion of the Earth’s center €(t) = r(t)/c in equation 174 about the midpoint to 
of the coherence time span Tyon: 


Co ¢(n) _ n 
E(t) = Eto) +) EE to" (177) 


Defining tg = 0 for convenience, the Earth’s orbital motion (excluding daily 
rotation) contribution to signal phase can then be written: 


ae i 
orp(t) = 2nKorn(t) «A (>. a | (178) 
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where m’ = min(m,N). 
It is also convenient to define new sky coordinates that capture the vector 
difference in signal phase (radians) between the source direction (a, 6) and 
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the detector’s direction from the Earth’s center (ap(to), dp) at time to: 


Nz (to) = 27 f(to)Te cos(d) cos(dp) cos[a — ap(to)], (181) 
Ny(to) = 27 f (to) Te cos(d) cos(dp) sin[a — ap(to)], (182) 
where Tg = Rp /c is the light travel time from the Earth’s center to the detec- 
tor. (See (Jaranowski and Krolak, 1999) for a similar sky coordinate definition.) 


Putting together equations 173-174 and 177-182 and absorbing phase con- 
stants into a single term $9, one obtains: 


eee k+1 
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Teoh 
+ me cos(2t) + ny (to) sin(2t), (183) 


where v‘“) are new coordinates, serving the role of effective frequencies and 
effective frequency derivatives: 
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where the insertion of powers of Toh is to make the coordinates dimensionless. 
Since all-sky searches to date have used only up to 1st-order frequency deriva- 
tives in first-stage analyis, it is useful to express v(tg) and v(t) = v (to) 
explicitly: 


Uta) = 2n=2" [flto) + Flto)E(to) A+ FElt0) A] (185) 


ito) = 2m ( Fen Ee OO) 45) «+ f8lt0)- A] (486) 

The form of equation 183 indicates the phase is linear with respect to the 
coordinates v“), nz and ny, which permits an analytic evaluation of the metric 
components (Pletsch, 2010) for a coherent search. Expressions appropriate for 
searching over a 2nd-order frequency derivative can be found in (Pletsch, 2010) 

Further, in the context of a semi-coherent search constructed from N.on 
coherently analyzed segments, one can systematically apply a refined metric 
in summing F-statistic values over the segments. In practice, a “coarse grid” 
for each segment 7 is defined by evaluating equation 158 to obtain the ooh In 
summing the F-statistic values, one must use a “fine grid” to avoid needless 
loss of SNR from signal evolution over the full observation period. As shown 
n (Pletsch, 2010), for the global correlation parameters, one can obtain the 
following approximation to the fine-grid metric from 


2 1 Uv] 
Jop = Neoh » Jap: (187) 
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Explicit evaluation of gag over many sidereal days leads to a fine grid that 
scales as ee for only the / coordinate (Pletsch, 2010), which is unsurprising, 
since the frequency derivative is the parameter driving the evolution of the 
frequency over time. Explicit expressions for gill and gag may be found in 
(Pletsch, 2010). One criticism (Wette, 2015) of this fine-grid metric approxi- 
mation, however, is that it does not explicitly take into account the changes in 
reference time implicit in each Taylor expansion for each segment. Nonethe- 
less, one finds empirically (Wette, 2016) that for semi-coherent searches the 
effective F-statistic mismatch grows much more slowly than implied by the 
Taylor expansion in Eqn. 177, allowing Eqn. 187 to be used successfully in 


Einstein@Home searches with large nominal metric mismatches. 


The Weave software infrastructure provides a more systematic approach to 
covering the parameter space volume in a templated search to ensure accept- 
able loss of SNR for true signals lying between template points (Wette ct al., 
2018). The Weave program combines together recent developments in tem- 
plate placement to use an optimal parameter-space metric (Wette and Prix, 
2013; Wette, 2015) and optimal template lattices (Wette, 2014). The package 
is versatile enough to be used in all-sky searches for unknown sources and 
in directed searches for particular sources, such as the Cas A and Vela Jr. 
supernova remnants (Abbott et al., 2022f). 


In brief, Weave creates a template grid in the parameter space for each time 
segment, a grid that is appropriate to computing the F-statistic for a coherence 
time Teon equal to the total observation period Tops divided by Necg. The 
spacing of the grid points in parameter space is set according to a metric (Wette 
and Prix, 2013; Wette, 2015) that ensures a worst-case maximum mismatch 
Meon defined by the fractional loss in summed F-statistic value due to a true 
signal not coinciding with a search template. 


Separately, a much finer grid is defined for the full observation period with 
respect to the midpoint of the observation period, one with its own mismatch 
parameter Msemi—coh, analogous to mon, Where the semi-coherent metric is 
the average of all the coherent metrics, which (unlike in the GCT approxi- 
mation) use a common reference time. The choice of the msemi—coh value is 
set empirically in a tradeoff between sensitivity and computational cost. The 
Weave package creates at initialization a mapping between each point in the 
semi-coherent template grid and a nearest corresponding point in each of the 
separate, coarser segment grids, accounting for frequency evolution. 


Finally, the discussion above has implicitly assumed analysis of data from 
a single detector. One may wonder if detection statistics based on two or more 
detectors require a finer template spacing, given the potential for better dis- 
crimination of signals by requiring coherent signal phase consistency among 
the detectors. For short coherence times there is indeed a finer discrimination 
from coherent summing when phase consistency is enforced and hence a need 
for finer sampling of frequency and sky location (Goetz and Riles, 2016). For 
much longer coherence times, however, this statement no longer holds. For 
example, the multi-detector F-statistic (Cutler and Schutz, 2005) has a co- 
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herent parameter space metric that is essentially unchanged from that of a 
single-detector F-statistic (Prix, 2007a). 

This perhaps surprising result can be understood from considering the in- 
trinsic motions of the detectors on the surface of an Earth in orbit. In order to 
maintain phase coherence for a single detector over the course of one day, one 
must track the detector’s relative motion around the Earth’s center a distance 
of order the diameter of the Earth (~13,000 km), larger than any detector 
pair separation. In addition, the Earth’s center travels a distance in its orbit 
of about 2.6 million km in one day, and more important to template spacing, 
deviates from a straight line by approximately 22,000 km. Given the phase fi- 
delity needed to account for these Earth-induced motions over coherence times 
much longer than this, the incorporation of additional detectors on the face 
of the Earth does not impose an extra burden on template placement. Note, 
though, that combining data coherently from Nget detectors of comparable 
sensitivity does improves SNR by the nominal desired Nae (Prix, 2007a). 


8.6.3 Viterbi methods and machine learning 


All of the search methods described so far use signal templates, explicitly or 
implicitly via favored frequency evolution. When searching a large parameter 
space volume with fine resolution, computational cost becomes formidable and 
often determinative of achievable sensitivity. Alternative approaches receiving 
increased attention rely upon more generic pattern recognition. 

The generic approach that has received most attention in recent years is 
based on Viterbi dynamical programming (Viterbi, 1967). To illustrate with 
a simplified example, consider finding a signal “trajectory” in a spectrogram, 
such as shown in Figure 9. A templated search might sum up the power for 
every possible trajectory allowed by the signal model and declare one or more 
candidate outliers based on a summed power of spectrogram pixels exceeding a 
pre-determined threshold. The Viterbi method (in its simplest form) dispenses 
with templates, seeking instead for the loudest trajectory that “moves” in time 
from left to right, where the degree of contiguity from one vertical column to 
the next is tunable. For example, a trajectory traveling from a pixel in column 
n and row m, to column n+1 may be constrained to change by no more than 
one row: |Mn41—M,y| < 1. For a trajectory that begins in row m, in column 1 
and travels to row my in the last column (NV), the number of possible trajec- 
tories is 3N—1. Maximizing the power over all possible such trajectories does 
not, however, require explicitly evaluating each power. The Viterbi algorithm 
leads to the insight that the trajectory with the highest summed power (for a 
strong enough signal) is also locally maximum, which allows rapid elimination 
of the vast majority of non-optimum trajectory segments and a remarkably 
fast evaluation of the detection statistic. 

The Viterbi method was first demonstrated in CW searches via a “spec- 
trogram” with each pixel representing a Bessel-weighted F-statistic evaluated 
over a 10-day period for Scorpius X-1 (Suvorova ect al., 2016) over the course 
of the initial LIGO S6 run (part of a Sco X-1 mock data challenge (Messenger 
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et al., 2015)). Follow-up analyses with additional refinements have been ap- 
plied to (Suvorova et al., 2017) or proposed (Melatos et al., 2021) for searches 
from the Advanced LIGO O1, O2 and O3 data (Abbott et al., 2017n, 2019e; 
Sun et al., 2020) (see section 4.3). Simultaneous tracking of stellar rotational 
phase and orbital phase (Melatos et al., 2021) offers a significant improvement 
in strain sensitivity relative to tracking of orbital phase alone (Suvorova ct al., 
2017). In addition, the Viterbi method has also been applied to searches for 
accreting millisecond pulsars (Middleton et al., 2020; Abbott et al., 2022d), 
isolated neutron stars (Sun et al., 2018; Millhouse et al., 2020; Abbott et al., 
2021) and for a post-merger remnant from the BNS merger GW170817 (Ab- 
bott et al., 2019d). The Viterbi method may see its largest gain in computation 
cost, though, from application to all-sky searches (Bayley et al., 2020; Abbott 
et al., 2022a). 


Although the hidden Markov Viterbi method has dramatic potential for 
reducing computational cost, it also has another important virtue; it is robust 
with respect to unknown and potentially stochastic frequency evolution that 
deviates from templated models. That flexibility makes the methodology es- 
pecially important for accreting systems like LMXBs (see section 2.1.4) and 
for extremely young sources, such as newborn neutron stars and post-merger 
hypermassive neutron stars (see section 3.10). 


Machine learning techniques, such as convolutional neutral networks, have 
received less attention, but offer similar gains in computational cost. One trains 
an algorithm on noise samples and signal++noise samples, for which machine 
learning detects an underlying pattern, producing an opaque but potentially 
effective algorithm for quickly yielding high detection statistic values for true 
signals. An early study (Dreissigacker et al., 2019) of single-detector data con- 
firms the enormous gain in computing cost possible, but does not suggest such 
automated algorithms achieve greater sensitivity. A follow-up study (Dreis- 
sigacker and Prix, 2020) examined machine learning on multi-detector data 
sets with realistic data gaps and non-Gaussian noise. Another study (Be- 
heshtipour and Papa, 2020) found that a convolutional neural network proved 
efficient in clustering Einstein@Home search outliers, to reduce computational 
cost in follow-up, with a different tuning found effective for identifying weak 
signals (Beheshtipour and Papa, 2021). Another recent study examined the 
potential for combining convolutional neural network analysis with Doppler 
demodulation for the Earth’s diurnal rotation in an all-sky search (Yamamoto 
and Tanaka, 2021). 


These generic methods are powerful in yielding rapid results, but require 
some care in use. For example, when searching a narrow band with instrumen- 
tal artifacts, the Viterbi method may seize upon the artifact and miss a nearby 
signal, although imposing consistency between different detectors can mitigate 
this problem (Bayley et al., 2020). An area of active research is understanding 
better the statistics of the loudest outlier in a Viterbi search, specifically, to 
understand the effective trials factor, a large value of which degrades strain 
sensitivity. In the event of a first detection via non-templated methods, there 
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remain, of course, fully templated methods available to assess more quantita- 
tively a candidate signal’s credibility and to estimate source parameters. 


3.7 Coping with non-Gaussian instrumental artifacts 


Non-Gaussian instrumental artifacts, especially spectral line artifacts, degrade 
CW searches. The degradation depends on the nature of the search. Stationary, 
narrow line artifacts generally do not significantly degrade targeted searches 
for known pulsars, for which long observation times permit extremely fine 
frequency resolution and known ephemerides permit that resolution to be ex- 
ploited. Periods during which a frequency-modulated signal overlaps with a 
known artifact can be vetoed or deweighted. On the other hand, an all-sky 
search is prone to contamination, especially in short data runs for which fre- 
quency modulation from certain sky directions may be limited, making a sta- 
tionary instrumental line resemble a signal template, at least in the first stage 
of a hierarchical search. 

For low assumed source spin-down (and no binary source modulation), the 
templates most prone to contamination lie near the ecliptic poles, where signal 
frequency modulation due to the Earth’s orbital motion would be small. At 
larger spin-down magnitudes, a stationary line can also lead to contamination 
of signal templates for which the frequency shift due to the Earth’s average 
acceleration toward the Sun largely cancels the assumed source spin-down. 
The associated templates tend to lie in a circular band concentric with the 
Sun’s average direction during the run with a radius and skyband thickness 
depending on the assumed frequency, spin-down and on the coherence time of 
the search (Abbott et al., 2008a). Such contamination is most pronounced for 
data runs short relative to a year. 

In principle, even a stationary line near an ecliptic pole should not be mis- 
taken for a true signal once a fully coherent algorithm is appled to assess that 
discrimination. The chance of an instrumental line displaying the residual fre- 
quency modulation (including that due to the Earth’s daily rotation) and asso- 
ciated phase modulation of a true signal is quite small. Moreover, the chance 
that two different detectors would display the same line artifact with pre- 
cisely the right time-dependent offset in phase to account for the daily change 
in relative positions of the detectors is quite small. For example, one veto 
method (Zhu et al., 2017) is based on turning off demodulation in the vicin- 
ity of an outlier template to determine if an even louder candidate is found. 
Another veto method, specific to the Frequency Hough search pipeline (see 
section 3.5.3), exploits characteristic patterns in the detection statistic varia- 
tion across search template parameter space created by stationary lines (Intini 
et al., 2020a). Similar considerations can be applied to following up outliers 
from Viterbi-based searches (Jones et al., 2022) (see section 3.6.3). Nonethe- 
less, lines are a major problem in CW searches because at initial stages of 
hierarchical searches, such discrimination is not available with tractable com- 
putational cost. Strong lines can trigger apparent loud signal outliers over 
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regions of parameter space, making outlier follow-up challenging. Simply ve- 
toing such a region because of a known contamination risks overlooking a true 
signal that would be recoverable in a deep search. 

Several methods have been developed for coping with these line-induced 
problems in early search stages (Leaci, 2015; Tenorio et al., 2021b), to reduce 
the burden of needless outlier followup while maintaining satisfactory detec- 
tion efficiency for true signals. The simplest method is to veto outliers known 
to be contaminated by a known line. This approach is effective in reducing 
computational cost, but does risk throwing away real and detectable signals. 
A more refined approach, one that need not rely upon prior knowledge of par- 
ticular lines is imposing consistency in signal strength seen in two or more 
detectors. For example, for two detectors of similar sensitivity one can require 
that individual detection statistic strengths in both detectors exceed a thresh- 
old and that the combined detection strength exceed both individual-detector 
strengths. Similarly, in a Bayesian approach one can impose consistency in the 
definition of the combined detections statistic (Keitel et al., 2014; Keitel and 
Prix, 2015; Keitel, 2016). An empirical background estimation to account for 
non-Gaussian contribution can be obtained (Isi ct al., 2020) via “sky-shifting,” 
that is, by evaluating template recovery strengths for identical source param- 
eters except for offsets in sky location. One can also require consistency in 
SNR across different data subsets for a putative outlier template, such as via 
a x? test for the separate contributions to the detection statistic (Sancho de la 
Jordana and Sintes, 2008). 

Another approach is “cleaning” of data prior to searching for CW sig- 
nal templates. Time-domain data cleaning has been used for general-purpose 
analysis (Allen et al., 1999; Meadors et al., 2014; Tiwari et al., 2015; Drig- 
gers et al., 2019; Davis et al., 2019; Vajente et al., 2020; Davis et al., 2021; 
Viets and Wade, 2021) where an auxiliary witness channel permits regres- 
sion of known noise. Such cleaning can remove both broadband and narrow 
contamination (Driggers ct al., 2019). A more CW-specific procedure can be 
carried out in the frequency domain in the absence of a witness channel — if 
a non-astrophysical source is clear. After creating DFTs one can replace bins 
known to be contaminated with randomly generated DFT coefficients consis- 
tent in magnitude with noise in neighboring bins (Abbott et al., 2009c). This 
approach potentially renders particular true signals less detectable or unde- 
tectable, particularly for sky locations near the ecliptic poles; hence injection 
simulations are needed to assess efficiency loss when setting upper limits in 
the abasence of a signal. 

Many spectral lines in a detector’s gravitational wave strain channel can 
be identified via correlation / coherence with lines observed in auxiliary chan- 
nels, such as for magnetometers or accelerometers, that monitor the environ- 
ment and that have no sensitivity to true astrophysical systems (Aasi et al., 
2015b; Covas et al., 2018). Others may not have a reliable witness channel, 
but come in “combs” of many lines with equal frequency spacings between 
adjacent lines, inconsistent with a plausible astrophysical source, allowing safe 
veto or cleaning (Goetz et al., 2021). Efficient tracking of known lines is an 
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active area of investigation, including tracking of lines that wander slightly in 
frequency (Daw et al., 2022). 

Traditionally, transient instrumental glitches in LIGO data that create nui- 
sances (sometimes severe) in searches for transient gravitational wave signal 
have not troubled CW searches much because their effect on overall noise 
level integrated over long time periods has been small. In the LIGO O3 data, 
however, a new class of extremely loud glitches with spectra peaking at low fre- 
quencies but visible as high as ~500 Hz appeared. These glitches of uncertain 
origin plagued both LIGO interferometers and occurred loudly and frequently 
enough to degrade sensitivity to CW signals in the low-frequency band. To 
cope with this new artifact, an ad hoc “self-gating” algorithm (Zweizig and 
Riles, 2021) was developed to taper the data in the time domain to zero dur- 
ing the affected intevals of ~seconds before creating DFTs for Fourier analy- 
sis. A more sophisticated, adaptive self-gating method (Steltner et al., 2022) 
achieved transient suppression with reduced deadtime. An earlier gating algo- 
rithm (Astone et al., 2005) was developed to cope with loud gltiches in initial 
Virgo data. 


3.8 Sensitivity depth 


A rough rule of thumb is convenient when assessing the detectability of a 
prospective CW signal for a given data set. Such a figure of merit is the 
sensitivity depth (Behnke et al., 2015). Its use arose in part because of the 
large variations in 1) methodologies with cost / sensitivity dependence on 
parameter space volume searched; 2) durations T.p, of data runs (or subsets) 
used in analyses; and 3) intrinsic detector sensitivity vs frequency. In part 
too, it avoids sometimes unwarranted assumptions based on idealized scaling 
with observation time. For example, a semi-coherent search with sensitivity 
improvements proportional to T a may require more computational resources 
than are available if Typ; becomes too large, especially since increasing Tops 
usually requires stepping more finely in parameter space. 

Instead, the sensitivity depth (Dreissigacker et al., 2018) addresses the 
“bottom line” with respect to a given intrinsic detector strain amplitude spec- 
tral noise density (square root of power spectral noise density S),): 


C=, (188) 


where ho is the quantity of interest, typically the 90% or 95% upper limit on a 
strain amplitude. By design the depth does not include a parametrized scaling 
with observation time. Hence the values for a given algorithm do depend on the 
particular data set. Reference (Dreissigacker et al., 2018) examines in detail 
the sensitivity depths achieved in searches of LIGO and Virgo data from the 
early initial LIGO $2 run to the first Advanced LIGO run O1. Values range 
for templated searches from ~1000 for targeted searches of ~2 years down to 
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~20 for the most sensitive all-sky search for CW signals in unknown binary 
systems. 


3.9 Upper limits and Sensitivities 


The CW search literature is rife with different conventions on how negative 
results (non-discoveries) are reported. This section gives a brief guide to the 
reader in understanding those variations and the reasons for them. 

Most analyses have produced frequentist upper limits at 95% (or 90%) 
confidence level, meaning that in a hypothetical ensemble of repeated exper- 
iments with the same underlying random noise contributions (but the same, 
non-random instrumental artifacts), a signal at the nominal upper limit value 
would have yielded a higher detection statistic 95% (90%) of the time. These 
upper limits are derived from or at least validated by simulated signals (injec- 
tions) and are quoted over narrow bands in frequency (usually 1 Hz or less), 
where wider bands necessarily have somewhat higher upper limits than most 
of the narrower bands from which they are composed. 

Deriving rigorous upper limits with extensive simulations in each individ- 
ual band is computationally expensive (particularly for 95% C.L.), so it has 
become common in recent years to derive instead “sensitivities” at, say, 95% ef- 
ficiency after following up and ruling out every outlier in each search band that 
lies above a nominal threshold (where the choice of threshold depends on a tar- 
get false alarm probability that varies considerably across different searches). 
These sensitivities are calibrated by deriving upper limits in a sparse sampling 
of narrow bands over the full search spectrum and finding an empirical scale 
factor between upper limits and average strain amplitude spectral densities 
for the data set, using a weighted average appropriate to the search. These 
sensitivities are not rigorous upper limits, particularly in disturbed bands, but 
give a useful interpretation of a non-detection. 

In highly disturbed bands with one or more strong instrumental lines, it is 
sometimes impractical to derive rigorous upper limits for some search methods 
or even to derive useful sensitivities. Such bands are vetoed and no upper limit 
quoted. As noted in section 3.7, when SFT cleaning of instrumental lines is 
used, one must take into account the resulting loss in detection efficiency in 
quoting upper limits. When strong lines are not vetoed or cleaned, upper limits 
in affected nearby bands typically suffer and may not apply at all to regions 
very near the ecliptic poles. 

Most quoted frequentist upper limits are population-averaged over the pa- 
rameter space searched, assuming random orientation of the stellar spin axis, 
and in the case of all-sky searches, random position on the sky. Detection ef- 
ficiency varies substantially for different angles of stellar inclination v (best 
efficiency for | cos(v)| near one, corresponding to circular polarization), and to 
a lesser extent over different regions of the sky. Because of this variation in 
sensitivity, the PowerFlux pipeline (see section 3.5.2) derives separate upper 
limits for circular polarization and linear polarization, where in each case the 
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95% C.L. upper limits are strict in the sense that 95% coverage is maintained 
separately for every position on the sky. Approximate population-averaged up- 
per limits can then be derived from the strict circular-polarization limits via 
multiplying by an empirically determined scale factor (typically ~2.3). 

As described in section 3.4.4, an alternative Bayesian analysis technique 
has been applied to targeted searches for known pulsars. In that approach a 
95% credible Bayesian upper limit on strain amplitude is obtained, which is 
interpreted as the analyst’s confidence that the true amplitude of a signal lies 
below that value, given the observed data and (conservative) prior beliefs in the 
parameter values. Bayesian notions of prior expectation have also influenced 
the construction of frequentist detection statistics. 


3.10 Transient CW sources 


In recent years, and particularly since the discovery of the binary neutron star 
merger GW170817, attention has turned to signal models that deviate from 
the canonical CW source of near-constant amplitude and very low intrinsic 
frequency evolution. Searches for two distinct classes of “near-CW” signals 
have been developed, one for sources of stable intrinsic frequency, but of large 
amplitude variations, and one for sources of rapid spin-down and concomitant 
amplitude decrease. The primary target motivating the first type of search is 
a neutron star glitch, in which a sudden stellar deformation appears, such as 
a ruptured crust, causing a sudden increase in the strength of gravitational 
waves emitted at twice the spin frequency of the star (Prix et al., 2011; Yim 
and Jones, 2020). The resulting stellar spin-down would be modest, leading 
to only small relative changes in frequency during the time required for the 
deformation to heal. Hence the search methods differ from “standard” CW 
methods primarily in allowing for a time-dependent strength. 

The danger in using the standard methods on a “transient CW” signal is 
that the data used prior to the glitch tends to reduce the integrated SNR, as 
does amplitude decay. To avoid this problem, an F-statistic-based method seg- 
ments the data and look separately for signals within individual segments and 
coherently or semi-coherently across different combinations of segments (Prix 
et al., 2011; Keitel, 2016; Keitel et al., 2019; Keitel and Ashton, 2018; Abbott 
et al., 2021h; Modafferi et al., 2021). 

Another class of near-CW source is a post-merger remnant, in which two 
neutron stars form a hypermassive neutron star (2-3 solar masses). Although 
one naively expects such a star to collapse promptly into a black hole, rapid 
rotation (rigid-body or differential) can delay the collapse for certain equa- 
tions of state (Baiotti and Rezzolla, 2017; Piro et al., 2017; Ravi and Lasky, 
2014). In extreme equations of state, the collapse may be delayed until the 
star’s rotation frequency has decreased dramatically (Ravi and Lasky, 2014). 
Given the enormous initial quadrupole asymmetry as two neutron stars begin 
to merge, one might hope for a substantial residual asymmetry in the min- 
utes, hours or even days during which a post-merger remnant persists. That 
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asymmetry might well lead to a rapid spin-down, one for which the truncated 
Taylor expansion in Equation 40 is a poor approximation. 

A recent search in LIGO data for a post-GW170817 remnant (Abbott 
et al., 2019d) used instead a model (for sensitivity determination) in which the 
frequency has an evolution similar to that of equation 5, but with a different 
normalization convention: 

dQ 

dt 
where 92 is the angular frequency of rotation, n is the braking index and 
k is a positive real constant. This equation leads to an explicit form for 
foaw(t) (Lasky et al., 2017b; Sarin et al., 2018): 


fa a. (190) 


(a+6)"" 


where Tgp is a characteristic time scale for spin-down: 
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Since the amplitude depends on frequency for fixed ellipticity (see equa- 
tion 14), one expects the amplitude to decrease monotonically too: 
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In addition, the product eJ,, is likely to decrease as the post-merger remnant 
spins down. 

A more significant hurdle to detection than fidelity of the signal model, 
however, is the typical distance at which binary neutron star mergers occur. 
GW170817 lay approximately 40 Mpc away, several orders of magnitude far- 
ther than the neutron stars sought in our own galaxy. The necessary ellipticity 
to generate a detectable signal is hence enormous; at the same time, such an 
ellipticity ensures a rapid-enough spin-down that no appreciable SNR could 
be achieved through integration over the signal’s duration at current detector 
sensitivities. Based on the total number of definitive BNS detections (two) (Ab- 
bott et al., 2017k, 2020a, 2021¢) during the O1 through O3 data runs and on 
the volume-time sampled in those runs, it appears that GW170817 was closer 
than the bulk of the BNS mergers expected in future runs. Detecting a CW 
signal from a post-merger remnant may require significantly more sensitive 
detectors than those that detected GW170817. Applicable search methods for 
such a rapidly evolving signal have been developed both well before (Thrane 
et al., 2011) and especially after (Thrane et al., 2011; Miller et al., 2018; Sun 
and Melatos, 2019; Oliver et al., 2019; Banagiri et al., 2019; Mytidis et al., 
2019; Miller et al., 2019a) the discovery of GW170817. 


ho 
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4 Results of continuous wave searches 


Searches have been carried out for continuous gravitational waves for five 
decades, starting with data from early detector prototypes (Levine and Steb- 
bins, 1972; Hirakawa et al., 1978; Livas, 1989; Suzuki, 1995). Although tran- 
sient gravitational wave discoveries to date have relied upon coincident signal 
detections in two or more detectors, a definitive continuous-wave source dis- 
covery can be accomplished, at least in principle, with a single gravitational 
wave detector. By definition, the source remains on, allowing follow-up veri- 
fication of the signal strength and of the distinctive Doppler modulations of 
signal frequency due to the Earth’s motion. In the event of an all-sky dis- 
covery, for which intrinsic sensitivity is necessarily limited by computational 
realities (see section 3.1), it is likely that a stable continuous signal could then 
a posteriori be detected in prior data sets via targeted searches. Hence a rela- 
tively large number of CW searches were carried out with both bar detectors 
and interferometer prototypes in the decades before the major 1st-generation 
interferometers began collecting data, as summarized in (Abbott et al., 2004). 


The most sensitive of the resulting early upper limits (Hirakawa et al., 1978; 
Suzuki, 1995; Astone et al., 2002b) came from bar detectors in their narrow 
bands of sensitivity. The Explorer detector reported (Astone et al., 2002b) an 
upper limit on spin-downless CW signals from the galactic center of 2.9 x 1074 
in a 0.06-Hz band near 921 Hz, based on 96 days of observation. A broader- 
band (~1 Hz) upper limit of 2.8x 10~2 was also reported (Astone et al., 2002a) 
from the Explorer detector based on a coherent 2-day search that allowing for 
stellar spin-down. In addition, searches for spin-downless CW waves from the 
galactic center and from the pulsar-rich globular cluster 47 Tucanae in two 
1-Hz bands near 900 Hz were carried out in Allegro detector data, yielding 
upper limits (Mauceli et al., 2000) of 8 x 10724. Finally, a narrowband (0.05 
Hz) search (Soida et al., 2003) was carried out with the TAMA interferometer 
near 935 Hz for continuous waves from the direction of Supernova 1987A, with 
an upper limit of 5 x 107? reported. 


When the initial LIGO interferometers and later the initial Virgo interfer- 
ometer began collecting data in the 2000’s, CW searches became more sensi- 
tive, both from improved detector sensitivity, and to a lesser extent, because 
search algorithms improved. In the following, brief summaries of the results 
from those searches will be given, with emphasis on results from searches in 
advanced detector data. As of this writing, numerous results from the third 
LIGO-Virgo observing run (O03) have appeared and will be featured where 
available, along with many results from the O1 and O2 runs, to illustrate the 
progression of sensitivities and algorithms during the Advanced LIGO and 
Virgo era to date. In recent years, research groups outside of the LIGO Sci- 
entific Collaboration, Virgo Collaboration and KAGRA Collaboration (LVK) 
have also carried out analyses of the public GW data, which is released ap- 
proximately 18 months after collection. Because of that delay, many additional 
results from the O3 data, beyond those described here, can be expected in the 
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coming months and years, perhaps in parallel with LVK results from the up- 
coming O4 data run (Abbott et al., 2020b). 


4.1 Targeted and narrowband searches for known pulsars 


In targeted searches for known pulsars using measured ephemerides from radio, 
optical, X-ray or y-ray observations valid over the gravitational wave obser- 
vation time, one can apply precise, well known corrections for phase of the 
signal, including modulations, because one knows the source phase evolution, 
its location and motion, the Earth’s location and motion, and the detector’s 
position and orientation on the Earth. 

Various approaches have been used in targeted searches in LIGO and Virgo 
data to date: 1) A time-domain heterodyne method (Dupuis and Woan, 2005) 
in which Bayesian posteriors are determined on the signal parameters that gov- 
ern absolute phase, amplitude and amplitude modulations (see section 3.4.4); 
2) a Fourier-domain determination of a “carrier” strength along with the 
strengths of two pairs of sidebands created by amplitude modulation from 
the Earth’s sidereal rotation of each detector’s antenna pattern (“5-Vector” 
method) (Astone et al., 2010b, 2012) (see section 3.4.5); and 3) a matched- 
filter method in which marginalization is carried out over unknown orienta- 
tion parameters (the “F-statistic”) (Jaranowski et al., 1998; Jaranowski and 
Krdélak, 2010) (see section 3.4.6). 

The first application of the heterodyne Bayesian method (Abbott et al., 
2004) in LIGO and GEO 600 S1 data (separately to each interferometer) led 
to upper limits on ho of a few times 10~?? for PSR J1939+2134 (for = 642 
Hz). Comparable upper limits were obtained from an implementation of the 
(frequentist) F-statistic (Abbott et al., 2004). Later applications of the het- 
erodyne Bayesian method incorporated a variety of improvements, including 
coherent treatment of multiple interferometers, marginalization over noise pa- 
rameters, a Markov Chain Monte Carlo search method for parameter estima- 
tion and joint searching over one and two times the stellar rotation frequency. 
At the same time the number of stars searched in each data run increased, 
along with closer partnership with radio and X-ray astronomers who provided 
ephemerides. In the $2 data, limits were placed on 28 pulsars, with a lowest 
strain limit of 1.7 x 10~74 (Abbott et al., 2005b). In the $3 and S4 data (an- 
alyzed jointly), limits were placed on 78 pulsars, with a lowest strain limit of 
2.6 x 10-7 (Abbott et al., 2007b). In the $5 data, limits were placed on 116 
pulsars, with a lowest strain limit of 2.3 x 10~?° (PSR J1603—7202) (Abbott 
et al., 2010). The lowest limit placed on ellipticity from the $5 search was 
7.0 x 10-8 (PSR J2124—3358). 

The final targeted-search results from initial LIGO and Virgo presented 
joint results from the LIGO $5 and S6 runs, and for the two low-frequency Crab 
and Vela pulsars, results from the Virgo VSR2 and VSR4 runs (Aasi et al., 
2014b). This synoptic paper presented results for 195 pulsars in total, where 
the lowest obtained strain limit was only slightly better than obtained from 
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the S5 data alone: 2.1 x 10-76 (PSR. J1910—5959D), with a lowest ellipticity 
upper limit of 6.7x 1078 (PSR J2124—3358). The use of Virgo VSR2 and VSR4 
data in this last analysis did, however, open up a new low-frequency spectrum, 
giving sensitivities approaching the spin-down limits for several pulsars other 
than the Crab, most notably the Vela pulsar, for which the spin-down limit 
was beaten (Abadie et al., 2011; Aasi et al., 2014b). The $5, S6, VSR2 and 
VSR4 analyses also included searches using the F-statistic and 5-vector algo- 
rithms applied to “high value” isolated pulsars for which the spin-down limits 
were approached or beaten. As expected, sensitivities obtained were compara- 
ble to those found in the Bayesian analysis. All three methods typically obtain 
somewhat better sensitivities when exploiting the inclination and polarization 
angles . and w inferred from pulsar wind nebulae observations for known pul- 
sars, such as Crab and Vela (for example, the F-statistic is refined to a more 
specific G-statistic (Jaranowski and Krolak, 2010)), although an unfavorable 
orientation can also lead to worse hg sensitivity. 


When Advanced LIGO data collection began in fall 2015 there was a sig- 
nificant improvement in broadband sensitivity and a dramatic improvement at 
the lowest frequency, thanks to improved seismic isolation (Aasi ct al., 2015a; 
Abbott et al., 2016c). The low-frequency improvements were, of course, helpful 
to the first binary black hole merger detection (Abbott et al., 2016b), but they 
also made a large number of known young pulsars accessible with respect to 
spin-down limit (energy conservation). Targeted searches were carried out in 
the O1 data using each search program (Abbott et al., 2017f), where method 
1) was applied to 200 stars, and methods 2) and 3) were applied to 11 and 10 
stars, respectively, for which the spin-down limit (equation 20) was likely to 
be beaten or approached, given the detector sensitivity. Results are shown in 
Figure 17, along with those from initial LIGO and Virgo searches. Highlights 
of these O1 searches included setting a lowest upper limit on strain amplitude 
of 1.6 x 10-76 (PSR. J1918—0642), setting a lowest upper limit on elliptic- 
ity of 1.3 x 10-8 (PSR J0636+5129) and beating the spin-down limit on 8 
stars (PSR J0205+6449, J0534+2200, J0835—4510, J1302—6350, J1813—1246, 
J1952+3252, J2043+2740, J2229+6114). Perhaps the most notable result was 
setting an upper limit on the Crab pulsar’s (PSR J0534+2200) energy loss to 
gravitational radiation at a level of 0.2% of the star’s total rotational enegy 
loss inferred from measured rotational spin-down. 


Similar searches were carried out for 221 known pulsars in the O1 and/or 
O2 data, with results summarized in Figure 17 (Abbott et al., 2019b). High- 
lights included beating the spin-down limit on 20 pulsars, a lowest upper limit 
on strain of 8.9 x 107-2” (PSR J1623—2631), a lowest upper limit on ellipticity 
of 5.8 x 10~° (PSR J0636+5129) and an upper limit on the Crab pulsar’s frac- 
tional energy loss to gravitational radiation of 0.02%. In addition, the upper 
limit on strain amplitude (1.5 x 107° for the MSP PSR J0711—6830 (frot = 
182 Hz) was only 30% above the star’s spin-down limit, corresponding to an 
ellipticity upper limit of 1.2 x 10-8. Upper limits are also presented in (Ab- 
bott et al., 2019b) on signals at the stellar rotation frequencies, along with 
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upper limits on the mass quadrupole moment Qa2 = €l,,4/ a (Ushomirsky 


et al., 2000). Initial analysis of the first six months of the LIGO and Virgo 03 
data set reduced further the upper limits on the Crab and Vela pulsar, along 
with those of three recycled pulsars, for which the spin-down limit has now 
been beaten (Abbott et al., 2020c). Similarly, a targeted search for the young, 
highly energetic star PSR, J0537—6910 (Abbott et al., 2021e) dived below the 
spin-down limit to an upper limit (95% CL) of 1 x 1076 for a GW frequency 
(123.8 Hz) of twice the rotation frequency. 

A separate analysis (Nieder et al., 2019) of the Advanced LIGO O1 and O02 
data for a newly discovered gamma-ray pulsar (PSR J0952—0607) also set an 
upper limit on emission amplitude of 6.6 x 10~?°. Upper limits on GW emission 
of amplitude 3.0 x 10~76 were also set on the black widow y-ray pulsar PSR 
J1653—0158 (Nieder et al., 2020) discovered in an Einstein@Home search. 

Recent cumulative results from targeted searches from the O1, O2 and O3 
data runs (Abbott et al., 20211) for 236 known pulsars in total are shown 
in Figure 19. Highlights include beating the spin-down limit on 23 pulsars, a 
lowest upper limit on strain of 4.7 x 107?” (PSR J1745—0952), a lowest upper 
limit on ellipticity of 5.26 x 10~9 (PSR J0711—6830), and an upper limit on 
the Crab pulsar’s fractional energy loss to gravitational radiation of 0.009%. 
In addition, the spin-down limit was beaten for two millisecond pulsars: PSR 
JO711—6830 (ho < 7.0 x 1072”, € < 5.3 x 107° for faw ~364 Hz) and PSR 
J0437—4715 (ho < 6.9x 1077", € < 8.5x10~° for faw ~347 Hz). These results 
also include a more general analysis searching simultaneously for a signal at 
one and two times the rotation frequency (Pitkin et al., 2015; Abbott et al, 
20211). Results from cumulative O1—O3a searches for seven additional pulsars 
were presented in (Ashok et al., 2021). 

The progressive improvement in noise level for the LIGO and Virgo de- 
tectors over the O1, O2 and O38 runs is reflected in Figures 17-19. Although 
more refined analyses have been brought to bear in parallel, the gains in astro- 
physical sensitivity come primarily from improving the instruments, for these 
targeted searches which already approach optimality. 

These upper limits assume the correctness of General Relativity in that 
antenna pattern calculations used in the searches assume two tensor polariza- 
tions in strain. Alternative theories of gravity can, in principle, support four 
additional polarizations (two scalar and two vector modes), which would lead 
to different antenna pattern sensitivities (Isi et al., 2015). Searches have been 
carried out for evidence of signals from the 200 targeted pulsars in the O1 
data exhibiting these other polarizations, using the heterodyned data prod- 
ucts. In no case was significant evidence of a non-standard signal seen, and 
upper limits were placed (Abbott et al., 2018a). 

The targeted-search upper limits in Figure 17 assume a fixed phase rela- 
tion between stellar rotation (measured by electromagnetic pulses) and gravi- 
tational wave emission (f; = frot). To allow for a more general scenario, such 
as slight differential rotation of EM- and GW-emitting regions, searches have 
also been carried out for signals very near in parameter space to those ex- 
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Fig. 17 Upper limits (95% CL) on ho for known pulsars from targeted searches in the LIGO 
O1 data (Abbott et al., 2017f) (closed stars). The gray band shows the a priori estimated 
sensitivity range of the search. Also plotted (closed squares) are the lowest upper limits 
from searches in initial LIGO and Virgo data and spin-down limits (closed triangles). Upper 
limits that lie below spin-down limits are outlined with a circle. 


pected from an ideal phase relation. These so-called “narrowband” searches 
allow a relative frequency deviation of ~10~. The first such search, using the 
F-statistic, was for the Crab pulsar in the Inital LIGO S$5 data set (Abbott 
et al., 2008b), which set a limit slightly below the Crab spin-down limit, de- 
spite a large trials factor of 3 x 10’, when the orientation of the assumed signal 
was aligned with observed Crab pulsar wind nebula X-ray jet axes (Ng and 
Romani, 2004), a limit five times higher than achieved in the same data set 
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Fig. 18 Upper limits (95% CL) on ho for 221 known pulsars from targeted searches in the 
LIGO O1 and/or O02 data (Abbott et al., 2019b) (closed stars). The pink band shows the 
a priori estimated sensitivity range of the search. Also plotted are spin-down limits (closed 
triangles). Upper limits that lie below spin-down limits are outlined with a circle. 


using a targeted search. A similar narrowband search was later carried out 
in initial Virgo VSR4 data for the Crab and Vela pulsars, using the 5-vector 
program, a search which yielded a Crab upper limit about two times below 
the spin-down limit and a Vela upper limit slightly higher than its spin-down 
limit (Aasi et al., 2015d). 

The 5-vector program was applied again to the Advanced LIGO O1 data 
set. Results from searches for 11 stars with expected sensitivities near the 
spin-down limits have been obtained from O1 data (Abbott et al., 2017e). 
In general, these limits are expected and found to be higher than the corre- 
sponding upper limits from targeted searches above because the increased pa- 
rameter space search implies an additional trials factor. Nonetheless, this first 
advanced detector narrowband search beat the spin-down limit on the Crab 
(PSR J0534+2200), Vela (PSR J0835—4510) and PSR J2229+6114. Later, a 
5-vector search of LIGO O2 data (Abbott et al., 2019c) for 33 known pulsars 
yielded the upper limits shown in Figure 20, along with the 11 (higher) O1 
upper limits. In this analysis, the spin-down limit was beaten for 8 known pul- 
sars, despite trials factors ranging from ~10° to ~10°. For the Crab pulsar, the 
strain upper limit was an order of magnitude lower than the spin-down limit, 
leading to a limit on fractional energy loss to gravitational waves of ~1%. 
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Fig. 19 Upper limits (95% CL) on ho for 237 known pulsars from targeted searches in 
the cumulative LIGO and Virgo O1—O3 data (Abbott et al., 2021l,c). The stars show 95% 
credible upper limits on the amplitudes of ho, while gray triangles represent the spin-down 
limits for each pulsar. canonical moment of inertia). For those pulsars which surpass their 
spin-down limits, their results are plotted within shaded circles. The pink curve gives an 
estimate of the expected strain sensitivity of all three detectors combined during the course 
of O3. The highlighted pulsars are those with the best ho, Q22 and spin-down ratio out of 
the pulsars which surpassed their spin-down limit, as well as the best ho limit out of the 
whole sample. 


Most recently, further sensitivity improvement was seen in O3 narrowband 
results (Abbott et al., 2021h), as shown in Fig. 21 for 18 known pulsars with 
spin-down limits within a factor of 3 of the expected sensitivity for which the 
spin-down limit is beaten for six pulsars. A separate analysis of O1—O38a data 
for seven other pulsars was carried out in (Ashok et al., 2021). 

Searches for accreting X-ray millisecond pulsars (AXMPs) (see section 2.1.4) 
require a modified narrowband approach in that nominal rotation frequencies 
are known, but with poor precision compared to that available for pulsars for 
which sustained monitoring is feasible. Their frequencies can vary rapidly (and 
likely with significant stochasticity) during active (accreting) phases and dur- 
ing quiescent phases can generally only be inferred. Given these uncertainties, 
including unknown stochastic contributions, a hidden Markov Viterbi method 
based on the 7-statistic!* has been applied to searches for five AXMPs in the 


14 The J-statistic is a weighted sum of powers from a large number of orbital sidebands 
generated by evaluating the F-statistic for a binary source, using a weighting governed by 
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02 LIGO data (Middleton et al., 2020) and to 20 AXMPs in the 03 LIGO 
data (Abbott et al., 2022d). The O3 search yielded strain amplitude sensitiv- 
ities in the range (5-24) x10~°, where estimated spin-down limits based on 
measured frequency derivatives lie in the range 10~78-10-?" with comparable 
to somewhat larger estimates based on torque balance (Abbott et al., 2022d) 
(see section 2.1.2). 


4.2 Directed searches for isolated stars 


Directed searches are those for which the source location is precisely known, 
but for which the signal’s gravitational wave phase evolution is unknown or 
poorly known. As discussed in section 3.6.1, the implied parameter space vol- 
ume of a truly broadband search will then depend sensitively upon the assumed 
age of the star. For a very young pulsar, one must search over not only the 
frequency and first frequency derivative (spin-down), but also over the second 
and possibly higher derivatives. 

Directed-search methods are also appropriate when searching for r-modes 
from known pulsars. The search band lies nominally near 4/3 the star’s rotation 
frequency, but has large systematic uncertainties of order 10% that depend on 
the unknown equation of state governing the modes (Idrisy et al., 2015; Caride 
et al., 2019). Hence, while the search band is much smaller than that for, say, 
a young neutron star with unknown frequency, the band is also much larger 
than that used in narrowband searches (see section 4.1), arguing for careful 
balancing of computational cost against sensitivity (Caride et al., 2019). 

The computational cost of fully coherent directed searches can be under- 
stood qualitatively from the scalings with coherence time implied by equa- 
tions 41-43, with more quantitative estimates based on the template place- 
ment considerations discussed in section 3.6.1. Semi-coherent searches have 
more complex scalings, but for long observation spans, generally achieve im- 
proved strain sensitivity with respect to fully coherent searches carried out 
over shorter subsets of the data set (which is typically necessary). 

The first such analysis in initial LIGO data used the F-statistic algo- 
rithm (Abadie et al., 2010) to search for the central compact object (X-ray 
source) at the center of the Cassiopeia A supernova remnant. As discussed 
in section 2, given the ~300-year presumed age of the star, one can derive a 
frequency-dependent upper limit on its strain emission of ~1.3 x 10774, as- 
suming its rotational energy loss has been dominated by gravitational wave 
emission. A coherent search was carried out in a 12-day period of LIGO S5 
data over the band 100-300 Hz, for which it was expected that the age-based 
limit could be tested with that data set (Wette ct al., 2008). The resulting 
upper limits did indeed beat the age-based limit over that band, reaching a 
minimum upper limit of 7 x 10775 at 150 Hz. That the limits were more than 
an order of magnitude higher than found in the full-S5 targeted searches for 


a set of Bessel functions Jn arising from the frequency modulation and incorporating the 
orbital phase of the binary system (Suvorova et al., 2017). 
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Fig. 20 Upper limits (95% CL) on ho for (11) 33 known pulsars from narrowband searches 
in the LIGO (O1) O2 data (Abbott et al., 2019c) (closed diamonds and triangles), where the 
GW frequency and derivative are allowed to vary by ~10~? with respect to the expectation 
from electromagnetic observations. For those pulsars known to have glitched in the O2 run, 
separate upper limits are shown for the epochs before the glitch (BG) and afterward (AG). 
Spin-down limits are shon as open circles, where error bars denote the uncertainties due 
to pulsar distances. Curves denote nominal sensitivities for the Ol and O2 runs for the 
individual LIGO Hanford (LHO) and Livingston (LLO) interferometers. 


known pulsars in that band reflected not only the much shorter observation 
time used (12 days vs. 23 months), but also the higher SNR threshold nec- 
essary to apply when searching over ~10!2 templates in f,, fz and f, for a 
300-year old star. 

This coherent approach over tractable intervals (Wette et al., 2008) was 
later applied to searches in the data from the last initial LIGO data run (S6) 
for nine young supernova remnants (Aasi et al., 2015e) and to a possible source 
at the core of the globular cluster NGC 6544 (Abbott et al., 2017m), achieving 
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Fig. 21 Upper limits (95% CL) on ho for known pulsars from narrowband searches in the 
LIGO O3 data (Abbott et al., 2021h). The red solid, blue dashed, and purple dotted curves 
show the expected sensitivities for H1, L1, and V1, respectively. The blue pentagons indicate 
the median 95% CL ULs from the 5n-vector search across all 10~* Hz sub-bands for each 
source. The black crosses indicate 95% CL ULs from the F-statistic search, which are set 
across the full search range for each target. The orange triangles indicate the spin-down 
limit, Aspin—down, With error bars that reflect uncertainty in the distance to each source. In 
a few cases the error bars are smaller than the size of the markers. 


upper limits on strain comparable to those found in the $5 data (lowest values 
ranging over ~(4-7)x10~?°, depending on source age (lower limits for older 
sources with lower trials factors from searching over frequency derivatives) 
The coherent F-statistic approach was applied to Advanced LIGO O1 
data (Abbott et al., 2019a) in a search for 15 supernova remnants and one 
nominal exoplanet with an unusual apparent orbit, which has been suggested 
to be a very nearby neutron star (Neuhduser et al., 2015) (see section 2.1.4). 
Best upper limits obtained ranged over ~(1-4) x 10~?°, depending on assumed 
source range. Figure 22 shows sample results for three of the supernova rem- 
nants, including Cas A, along with that for Fomalhaut b. In this analysis, 
separate “deep” and “wide” analyses were applied to three of the supernova 
remnants, including Vela Jr., to account for large uncertainties in source age, 
where deep searches could be carried out for older sources, requiring a smaller 
range in frequency derivatives. A similar approach was used to probe the O2 
data for 12 supernova remnants, restricting attention to frequencies below 150 
Hz, applying coherence times ranging from 12 to 55.9 days (Lindblom and 
Owen, 2020). A recent coherent search of 8.7-day and 12.8-day subsets O2 
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Fig. 22 Upper limits (95% CL) on ho (dots) for 3 supernova remnant cores and nominal 
exoplanet (but possible neutron star) Fomalhaut b, using coherent F-statistic searches of 
O1 data (Abbott et al., 2019a). Upper left: G1.9+0.3; upper right: Vela Jr.; lower left: 
Cassiopeia A; lower right: Fomalhaut b. The horizontal lines indicate nominal age-based 
limits (equation 28). 


LIGO data (Owen et al., 2022) for CW radiation from a Supernova 1987A 
remnant beat the age-based indirect limit. A Viterbi-based search for Fomal- 
haut b was applied to O02 data (Jones and Sun, 2021). 


As one might imagine, a semi-coherent approach has the potential to im- 
prove upon a single coherent directed search. One demonstration of the method 
in initial LIGO S5 data searched for a source at the galactic center (Aasi et al., 
2013a), using 630 segments of 11.5 hours each, where F-statistic values aver- 
aged over the segments were computed, where the global correlation transform 
template mapping was used in combining the F-statistic values over the seg- 
ments. A similar but more sensitive semi-coherent approach was applied in 
a computationally intensive Einstein@Home (see section 4.4) S6 search for a 
CW signal from Cas A (Zhu et al., 2016). This search used 44 segments of 140 
hours, again applying the global correlation transform template gridding and 
summing. 

The same method was applied to an Einstein@Home search in Advanced 
LIGO O1 data for three supernova remnants: Cas A, Vela Jr. and G347.3 (Ming 
et al., 2019). Figure 23 shows the results of the three searches, together with 
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results from the coherent search of a subset of the same data set (Abbott 
et al., 2019a). The semi-coherent search which exploits the full data set, typi- 
cally achieves a factor of two improvement in strain sensitivity over the coher- 
ent search over a data subset. This search also applied a search optimization 
method (Ming et al., 2016a) to choose coherence times and segmentations 
for each source, where the optimization attempts to take into account rela- 
tive probabilities for detection, given available astronomical information. In 
this instance the segmentations chosen were 12 245-hour segments (Cas A), 
8 369-hour segments (Vela Jr.) and 6 489-hour segments (G347.3). As seen 
in Figure 23, best upper limits obtained were ~1 x 107°. Another semi- 
coherent directed search for the galactic center, based on the Frequency Hough 
method (Antonucci et al., 2008), accelerated by the Band-Sampled Data (BSD) 
use of DFTs (Piccinni et al., 2018), was carried out using the Advanced LIGO 
O02 data (Piccinni et al., 2020). Another O2 analysis (Ming et al., 2022) 
searched for a signal from the G347.3 using a semi-coherent F-statistic im- 
plementation in Einstein@Home. 

Several distinct semi-coherent directed searches have been carried out to 
date using O3 data. First came results from three methods applied to 15 
supernova remnants using the O3a data (Abbott et al., 2021k), one method 
being a semi-coherent, BSD-accelerated Frequency Hough search and the other 
two methods being less sensitive but more robust Viterbi (Sun et al., 2018) 
searches. The two Viterbi methods searched for either signal at only a single 
frequency (Sun et al., 2018) (assumed to be twice the unknown rotation fre- 
quency) or at both one frequency and its doubled value (Sun et al., 2019). 
All three methods were applied to seven stars, with only the single-frequency 
Viterbi method applied to another eight stars. The results from the Frequency 
Hough search are shown in Fig. 24. 

A separate O3a analysis (Abbott et al., 2022f) used the Weave implementa- 
tion (see section 3.6.2) of a semi-coherent F-statistic search. While the package 
is versatile enough to be used in all-sky searches for unknown sources, a sim- 
pler configuration, applicable to well localized sources, was used to search in 
the O8a data for the Cas A and Vela Jr. supernova remnants. Figure 25 shows 
the results in comparison with earlier searches for these two sources in O1, O2 
and O8a data. 

These 95%-efficiency sensitivities to Cas A and Vela Jr. can be translated 
into sensitivities to ellipticity, as shown in Fig. 26. The quadratic dependence 
of strain on frequency for fixed ellipticity (see Eqn. 14) leads to dramatically 
better sensitivity to ellipticity at higher frequencies, reaching as low as € = 
2x 107-8 near 1000 Hz for the more optimistic assumption of Vela Jr. distance 
(0.2 kpc). 

Another approach (Dhurandhar et al., 2008) for directed searches is based 
on cross correlation of independent data streams. The most straightforward 
method defines bins in detector-frame frequency and uses short coherence 
times, as in directional searches for stochastic gravitational radiation, (Ballmer, 
2006a; Abbott et al., 2017b) which can be used to search for both isolated and 
binary sources, albeit with limited sensitivity. One can use finer frequency 
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Fig. 23 Upper limits (90% CL) on ho for Cassiopeia A, Vela Jr. and G347.3 (dots) 
from semi-coherent Einstein@Home directed F-statistic searches in Advanced LIGO O1 
data, (Ming et al., 2019; Papa et al., 2020), shown with previous (higher) coherent-search 
limits using subsets of the Ol data. The dashed curves denote estimated 95% CL upper 
limits based on the 90% CL values. 


binning, however, when correcting explicitly for Doppler modulation of the 
signal. Cross-correlation methods are especially robust against wrong assump- 
tions about phase evolution and are attractive in searching for a very young 
object, such as a hypothetical neutron star remaining from Supernova 1987A 
(see (Ashton et al., 2017) for a discussion of potential degradation of coherent 
searches from neutron star glitches, (Page et al., 2020) for evidence of a hidden 
star from an excess of infrared emission and (Greco et al., 2022) for evidence 
of pulsar wind nebula). A cross-correlation search for SN 1987A, including de- 
modulation for effects from the motion of the Earth, (Chung et al., 2011) was 
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Fig. 24 Sensitivity estimates (95% efficiency) 
Frequency Hough (BSD-accelerated) search (Abbott et al., 2021k). The dotted curves rep- 
resent the estimated ho in the full band of H, L and V detectors searched by the pipeline. 
The crosses represent the frequentist strain upper limits at 95% confidence level obtained 
empirically in sample sub-bands of 1 Hz. Horizontal lines are the indirect age-based limit 
(Eqn. 28). The limit is beaten across the full band also using Virgo data, except for the 
most disturbed regions, for G65.7+1.2, G189.1+3.0 and G266.2-1.2/Vela Jr. The remaining 
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Fig. 25 Top panel: Estimated gravitational wave strain amplitude sensitivities (95% effi- 
ciency) in each 0.1-Hz sub-band for the O3a Cas A (red band) and Vela Jr. (cyan band) 
searches (Abbott et al., 2022f). Conservative uncertainty bands of +7% are indicated, to 
account for statistical and systematic uncertainties in estimating sensitivity depths, includ- 
ing calibration uncertainties. Black triangles (upright — Cas A, inverted — Vela Jr.) denote 
0.1-Hz bands for which rigorous upper limits are used to determine estimated sensitivity vs. 
frequency. Additional results from prior searches for Cas A and Vela Jr. are also shown: O1 
Einstein@Home 90% C.L. upper limits for Cas A (magenta curve) and for Vela Jr. (green 
curve) (Ming et al., 2019); O3a Cas A and Vela Jr. 95% C.L. upper limits using a model- 
robust Viterbi method (orange curve) (Abbott ect al., 2021k); O3a Vela Jr. 95% C.L. upper 
limits using the template-based Frequency Hough method (black curve) (Abbott et al., 
2021k). The solid red horizontal line indicates the age-based upper limit on Cas A strain 
amplitude. The dashed (dotted) horizonal blue lines indicate the optimistic (pessimistic) 
age-based upper limit on Vela Jr. strain amplitude, assuming an age and distance of 700 
yr and 0.2 kpc (5100 yr and 1.0 kpc). Bottom panel: Magnification of the sensitivity bands 
from the O3a Weave search over most of the search band (~40-976 Hz), with 1-o statistical 


uncertainties shown for the individual sparsely sampled upper limits used to estimate the 
depth. 
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Fig. 26 Estimated ellipticity sensitivities (95% efficiency) in each 0.1-Hz sub-band for the 
O3a (Weave-based) Cas A (red) and Vela Jr. (blue, magenta) searches, derived from the 
strain amplitude sensitivities shown in Fig. 25 assuming a source distance of 3.3 kpc for Cas 
A, and assuming source distances of 1.0 kpc and 0.2 kpc for Vela Jr. (color online). 


carried out in initial LIGO data (Sun et al., 2016) Recent application of cross- 
correlation methods to directed searches for binary sources will be discussed 
in the next section (4.3). 


Directed searches for particular sources require making choices, that is, to 
prioritize among a wide set of potential targets in deciding how best to apply 
computational resources and analyst time. Recent work (Ming et al., 2016b, 
2018) has taken a probabilistic approach to address this problem, based on 
source age and distance information (including sometimes large uncertainties) 
along with detector sensitivity, an approach that may be generalized to pa- 
rameter choices in both directed and all-sky searches. 


Searches for r-modes radiation from known pulsars are less challenging 
computationally than truly broadband directed searches, because the range of 
expected frequencies is better known. Nonetheless there is substantial theoret- 
ical uncertainty in the ratio between GW emission frequency and rotation fre- 
quency. Although the nominal ratio is 4/3 in the slow-spinning, non-relativistic 
regime, there are substantial corrections for fast-spinning stars and for stellar 
compactness that depend on the equation of state (Yoshida et al., 2005; Jasi- 
ulek and Chirenti, 2017; Idrisy et al., 2015; Caride et al., 2019), leading to a 
significant range in possible ratios. Following (Caride et al., 2019), the ratio 
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can be written: 


2 
fow . 4 B( Prot ) ; (193) 


frot fikepler 


where fKepler is the Kepler frequency of the star (rotation frequency at which 
centrifugal forces destroy the star), A is a parameter dependent on the equation 
of state, with an estimated allowed range of 1.39-1.57 (Idrisy et al., 2015) and 
B is acorrection term for high spin with an estimated maximum value of Bmax 
= 0.195 (Caride et al., 2019). 

Using these assumptions, several searches have been carried out explicitly 
for such r-modes: 1) an analysis of Ol and O2 data for emission from the young, 
energetic pulsar PSR, J0537—6910 (Fesik and Papa, 2020) (see section 2.1.4, 
which reached to within an order of magnitude of the strain spin-down limit; 
2) an analysis of Ol and O2 data for emission from the younger, compara- 
bly energetic and much closer Crab pulsar, for which the spin-down limit was 
surpassed by an order of magnitude (Rajbhandari et al., 2021); and a recent 
search in the O03 LIGO and Virgo data (Abbott et al., 2021d) which placed 
stringent constraints on theoretical models for r-mode-driven spin-down in 
J0537—6910, especially for higher frequencies for which upper limits reach be- 
low the spin-down limit. These latter results which attempt to address directly 
the evidence for r-modes in inter-glitch JO537-6910 spin-down are shown for 
the frequency band of the search (86-97 Hz) in Fig. 27. 


4.3 Directed searches for binary stars 


For known binary pulsars with measured timing ephemerides, targeted searches 
work well, and upper limits have been reported for many stars, as described 
in section 4.1. But searching for known (possibly accreting) neutron stars in 
binary systems not exhibiting pulsations or for entirely unknown stars in bi- 
nary systems once again significantly increases the parameter space, relative 
to the corresponding isolated star searches, posing new algorithmic challenges 
and computing costs. 

Searches for Sco X-1 in O1 data were carried out with several methods: 
1) a “Sideband” method (Messenger and Woan, 2007; Sammut et al., 2014; 
Suvorova et al., 2016; Abbott et al., 2017n; Sun, 2018) based on summing 
power in orbital sideband frequencies; 2) a non-demodulated cross-correlation 
method (Ballmer, 2006b; Abbott et al., 2017b) and 3) a demodulated cross- 
correlation method (Whelan et al., 2015; Abbott et al., 2017p; Meadors et al., 
2018). The demodulated cross-correlation method has proven to be the most 
sensitive method to date in such searches on a fixed data set for templated 
signal models without stochasticity, as expected from a previous mock data 
challenge (Messenger et al., 2015) including these methods and others (Goetz 
and Riles, 2011, 2016; Meadors et al., 2016, 2017; van der Putten et al., 2010), 
and as shown in Figure 28. Computationally intensive methods using the F- 
statistic, however, may eventually improve upon it (Leaci and Prix, 2015). 
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Fig. 27 Upper limits on GW amplitude ho obtained from searches for r-modes emission 
from PSR J0537-6910 in the O3 LIGO data using the F-statistic/G-statistic and 5-vector 
methods (Abbott et al., 2021d). The shaded band indicates the full range of results of 
the F-statistic/G-statistic pipeline. The dashed lines are defined by the stiffest and softest 
equations of state considered in the analysis and enclose a range of theoretical ho. 


Follow-up Sco X-1 searches of the O2 data were based on the Viteri method 
using the J-statistic (Abbott et al., 2019e, 2022e). 


One complication in Sco X-1 searches is potential spin wandering due to 
fluctuations in accretion from its companion (Mukherjee et al., 2018), which 
limits the length of a coherence time that can be assumed safe for a signal 
template. One previous fully coherent search (Aasi et al., 2015c) restricted 
its coherence length to 10 days, to be conservative. Semi-coherent and cross- 
correlation methods (Suvorova et al., 2016; Meadors et al., 2016; Ballmer, 
2006a; Whelan et al., 2015) should be more robust against wandering. Fig- 
ure 29 shows results from the recent Viterbi-based Sco X-1 search (Abbott 
et al., 2022c) in the O3 data using the J-statistic (Suvorova et al., 2017), in 
which results for different assumptions about Sco X-1 orientation are made. 
The implied limits on intrinsic strain amplitude ho are lowest in the most fa- 
vorable case of circular polarization, less favorable for an inclination angle 
= 44° consistent with observations of its radio lobes (Fomalont et al., 2001), 
and least favorable for a strain amplitude marginalized over unknown incli- 
nation. Figure 30 shows a comparison of Sco X-1 upper limits (marginalized 
over the unknown stellar inclination angle) obtained from a variety of meth- 
ods applied to the O1, O2 and O03 LIGO data, with comparisons to the torque 
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Fig. 28 Upper limits (95% CL) on ho for Scorpius X-1 from Advanced LIGO O1 data, us- 
ing several different search methods: a “radiometer” search using stochastic analysis meth- 
ods (Abbott et al., 2017b) and fine frequency binning, a Viterbi method based on a Bessel- 
weighted F-statistic (Abbott et al., 2017n) and a templated cross-correlation method (Ab- 
bott et al., 2017p). The dashed line indicates the torque-balance benchmark defined in 
equation 33 for accretion at the stellar radius. 


balance limit for a stellar vs. Alfvén radius for accretion (see section 2.1.2). 
The O38 inclination-averaged strain upper limits (Abbott et al., 2022c) shown 
in Figure 30 do not quite reach as low as the torque-balance benchmark in 
equation (33), but when results from the computationally more demanding 
CrossCorr search of the 03 data are available, one may expect the benchmark 
to be reached at the lower frequencies (below ~300 Hz). As Advanced LIGO 
sensitivity continues to improve and with longer data runs, future searches 
should progressively probe to higher frequencies along this benchmark. 

Possessing more definitive information on the rotation frequency of Sco X-1 
could potentially make the difference between missing and detecting its gravi- 
tational waves in advanced detector data, by both permitting longer coherence- 
time searches and reducing the statistical trials factor and thereby the thresh- 
old needed to identify an interesting outlier. More intensive measurements 
and analysis of Sco X-1 X-ray emission could yield a dramatic scientific pay- 
off (Galaudage et al., 2021). 

Until recently, CW searches for known LXMB systems focused almost ex- 
clusively on Scorpius X-1, although (Meadors et al., 2017) did also include 
limits from narrowband searches around three particular frequencies of in- 
terest for XTE J1751-—305, given X-ray observations of a potential r-mode 
excitation(Strohmayer and Mahmoodifar, 2014a). More attention is turning 
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Fig. 29 Frequentist effective wave strain upper limits at 95% confidence from LIGO O3 
data (Abbott et al., 2022e) as a function of sub-band frequency, for three scenarios: circular 
polarization with . = 0 (blue stars), v & 44° based on radio observations (see (Fomalont 
et al., 2001); orange dots), and a flat prior on cose (green dots). Indirect torque-balance 
upper limits (see Section 2.1.2) for two torque lever arms are also shown: the stellar radius 
(red solid line) and the Alfvén radius (dashed red line). 


now to other accreting systems, such as Cygnus X-2 (Premachandra et al., 
2016; Galaudage et al., 2021). In addition, recent searches were carried out in 
Advanced LIGO O2 data for five systems (Middleton et al., 2020) and in 03 
data for 20 accreting millisecond pulsars (Abbott et al., 2022d), both using 
a Viterbi hidden Markov method (Suvorova et al., 2017) and both exploiting 
the relatively good precision with which the stellar rotation frequencies are 
known (see section 4.1). 


4.4 All-sky searches for isolated stars 
4.4.1 Overview of search pipelines in use 


Various semi-coherent algorithmic approaches have been tried, many based 
in some way on the “Stack Slide” algorithm (Brady et al., 1998; Brady and 
Creighton, 2000; Cutler et al., 2005; Mendell and Landry, 2005) in which the 
strain powers from Fourier transforms computed over each coherently ana- 
lyzed segment are stacked on each other after sliding each transform some 
number of bins to account for Doppler modulation of the source frequency 
(see section 3.5.1). One algorithm is a direct implementation of this idea called 
StackSlide (Mendell and Landry, 2005). 
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Fig. 30 Comparison of 95% CL upper limits on ho from Sco X-1 emission from multiple 
search methods carried out in O1, O2 and O3 data: black — Ol CrossCorr search (Abbott 
et al., 2017p), brown — O2 CrossCorr search (Zhang et al., 2021), pink - O1-O3 Radiometer 
search (Abbott et al., 2021j); and green — O03 Viterbi search (Abbott et al., 2022c). The 
indirect torque-balance upper limits (see Section V C), using the stellar radius (red solid 
line) and the Alfvén radius (dashed red line), are also plotted, marginalized over stellar 
inclination angle. 


Other implementations (Krishnan et al., 2004; Antonucci et al., 2008) are 
based on the Hough transform approach, (Hough, 1959, 1962) in which for each 
segment a detection statistic is compared to a threshold and given a value of 0 
or 1. The unity values were later refined to be adaptive non-unity weights, to 
account for variations in noise and detector antenna pattern (Palomba et al., 
2005; Sintes and Krishnan, 2006). The sums of those weights are accumu- 
lated in parameter space “maps,” with high counts warranting follow-up. The 
Hough approach offers greater computational efficiency from reducing floating 
point operations, along with robustness against non-Gaussian artifacts (Ab- 
bott et al., 2008a) (see section 3.7). The Hough approach has been imple- 
mented in two distinct search pipelines, the “Sky Hough” (Krishnan et al., 
2004; Abbott et al., 2005a) and “Frequency Hough” (Antonucci et al., 2008; 
Astone et al., 2014a; Aasi et al., 2016a) programs, named after the different 
parameter spaces chosen in which to accumulate weight sums. 


Another implementation, known as PowerF lux, (Abbott et al., 2008a; Der- 
gachev, 2005; Dergachev and Riles, 2005; Dergachev, 2010a,b, 2013) improves 
upon the StackSlide method by weighting segments by the inverse variance of 
the estimated (usually non-stationary) noise and by searching explicitly over 
different assumed polarizations while including the antenna pattern correction 
factors in the noise weighting (see section 3.5.2). 
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Yet another method uses coincidences among F-statistic outliers (see sec- 
tion 3.4.6) in multiple time segments typically longer than those used in the 
semi-coherent approaches (Astone et al., 2010a; Aasi et al., 2014c), where the 
implementation is carried out in the time domain (hereafter denoted as TD-F- 
statistic), with systematic follow-up of outliers carried out through progressive 
increase of coherence time (Sieniawska et al., 2019). 


The deepest wideband searches (including wide in frequency derivative 
range) achieved to date in given fixed data sets have stacked F-statistic val- 
ues over time segments semi-coherently (see section 3.5.4) and have used the 
resources of the distributed computing project Einstein@Home (Abbott et al., 
2009c) based on the same software infrastructure (BOINC) (University of Cal- 
ifornia, 2002) developed for the Seti@Home project (Anderson et al., 2002). 
Einstein@Home encourages volunteers to download narrow-band segments of 
LIGO data and carry out a semi-coherent F-statistic search over a small patch 
of sky. Results are automatically returned to an Einstein@Home server and 
recorded, with every set of templates analyzed independently by host com- 
puters owned by at least two different volunteers. Einstein@Home Scientists 
then carry out post-processing to follow up on promising outliers found. This 
project has been remarkably successful in engaging the public (hundreds of 
thousands of volunteers and 750,000 host computers to date) in forefront sci- 
ence while making good use of idle computer cycles to carry out searches that 
would otherwise exceed the capacity of dedicated gravitational wave comput- 
ing clusters. 


The availability of the Einstein@Home platform has driven the evolution 
of semi-coherent stacked F-statistic techniques. This evolution has led to im- 
creased search sophistication and sensitivity over the last decade and a half, in 
general, including for related pipelines outside of that distributed computing 
framework, such as Weave (which has a memory footprint incompatible with 
Einstein@Home). Particular improvements have included search setup opti- 
mization (Cutler et al., 2005; Prix and Shaltev, 2012; Shaltev, 2016), more ef- 
ficient semi-coherent stacking and template placement, (Prix, 2007a,b; Pletsch, 
2008; Pletsch and Allen, 2009; Wette and Prix, 2013; Wette, 2014, 2015, 2016; 
Wette et al., 2018; Walsh et al., 2019) automated vetoing of instrumental 
lines, (Keitel et al., 2014; Keitel and Prix, 2015; Keitel, 2016) and hierarchical 
outlier followup and veto (Shaltev et al., 2014; Papa et al., 2016; Singh et al., 
2017; Zhu et al., 2017; Ashton and Prix, 2018; Intini et al., 2020b). 


Technical challenges in distributed computing include efficient data transfer 
to/from host computers and running on many computing platforms of greatly 
varying CPU and memory capabilities. The large computing resources avail- 
able via distributed computing can be used to enlarge the parameter space 
searched or to probe more deeply in the noise than is feasible on current com- 
puting clusters, but optimization must account for scaling of computing cost 
with the target range of frequency and frequency derivative and weigh the 
benefit of longer coherence time for sensitivity against the incurred cost (see 
section 3.1). 
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The F-statistic-stacking techniques can also be used, of course, in less 
powerful computing environments, with different tunings, e.g., shorter coher- 
ence times per segment. These techniques can also be used for systematic 
follow-up of outliers found in first-stage semi-coherent F-statistic searches or 
in searches using other semi-coherent methods (Walsh et al., 2019), including 
both all-sky and directed searches. One general-purpose, multi-stage approach 
uses the python wrapper PyFstat for F-statistic summing (Ashton and Prix, 
2018; Keitel et al., 2021) and a Markov Chain Monte Carlo search through 
parameter space to “zero in” on signals (Tenorio ct al., 2021a). This method 
systematically lengthens segment coherence times (hence reducing segment 
counts per observational run) simultaneously with narrowing of the parameter 
space volume, while guided by the parameters of the loudest survivors from 
each stage. 

A comparison of many of these all-sky search methods was carried out via 
a mock data challenge using initial LIGO data, (Walsh ct al., 2016) and these 
methods have been applied to searches of the Advanced LIGO O1-O3 data 
sets (Abbott et al., 2017a,d, 2018b, 2019a; Steltner et al., 2021; Abbott et al., 
2021b, 2022a). Unsurprisingly, the all-sky search enabled by Einstein@Home 
computing resources displayed consistently better sensitivity than the other 
methods in the mock data challenge. 

A newcomer all-sky search pipeline, known as the SOAP pipeline (Bayley 
et al., 2020), uses a Viterbi approach to seek trajectories in spectrograms for 
which each time segment is represented by the average spectrum over a 24- 
hour period using a 30-minute coherence time. Although not as sensitive as 
the pipelines described above, the technique is blazingly fast, in comparison, 
offering the potential of rapid discovery for observing runs with much improved 
detector noise. Perhaps more important, because the algorithm is untemplated, 
it has the additional potential of detecting new (strong) signals that do not 
follow the models sought by other isolated-star pipelines, including long-period 
binary systems. 


4.4.2 Results from all-sky, isolated-star searches of LIGO and Virgo data 


The Sky Hough algorithm was used to produce all-sky upper limits in the 
200-400 Hz band of the LIGO S2 data (Abbott et al., 2005a), based on a total 
of 3800 30-minute segments of data from the three LIGO interferometers. 
The StackSlide, Sky Hough and PowerFlux methods were used to produce 
all-sky upper limits in the 50-1000 band of the LIGO S4 data (Abbott et al., 
2008a). The first Einstein@Home all-sky search was carried out too on the S4 
data (Abbott et al., 2009c). 

The PowerF lux algorithm was used to produce all-sky upper limits in the 
50-1100 Hz band of the first eight months of LIGO $5 data (Abbott et al., 
2009a). The sheer length of data for the full 23-month $5 run required sub- 
stantial upgrade of the program which was then used to produce all-sky upper 
limits in the 50-800 Hz band of the full data set, based on a total of more 
than 80,000 (50%-overlapped) 30-minute segments from the H1 and L1 data. 
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This PowerFlux result (Abadie et al., 2012) included a three-stage hierarchical 
search with a follow-up procedure of loud candidates based on loose coherence 
(see section 3.3.4). A Sky Hough search of the S5 data consisted of a coin- 
cidence analyis of data sets from two separate approximately 1-year subsets 
of the data over the 50-1000 Hz band. Einstein@Home too was applied in se- 
quential analyses to the early $5 (Abbott et al., 2009b) and to the full S5(Aasi 
et al., 2013b). A final all-sky initial LIGO PowerFlux analysis of the S6 data 
set (Abbott et al., 2016a) included a 5-stage hierachical search with longer and 
longer effective coherence times over 100-1500 Hz within the loose coherence 
framework. The S6 Einstein@Home search (Abbott et al., 2016) achieved the 
most sensitive all-sky results from any of the initial LIGO data sets, reaching 
upper limit values as low as 5.5 x 10729. 


When initial Virgo VSR1 data became available, a direct time-domain 
implementation of the F-statistic (Astone et al., 2010a) was applied to a search 
of it for the 100-1000 Hz band (Aasi et al., 2014c). Later, the Frequency Hough 
method was applied to data from the initial Virgo VSR2 and VSR4 runs over 
the 20-128 Hz band, the first time an all-sky search was applied to frequencies 
below 50 Hz (Aasi et al., 2016a). 


Since Advanced LIGO observing has begun, multiple all-sky search pro- 
grams have been applied to data from the first three observing runs, O1, O2 
and O3. The first publications based on O1 data focused on lower frequencies. 
Four pipelines (PowerF lux, Sky Hough, Frequency Hough and TD-F-statistic) 
covered the band 20-475 Hz and a spin-down range [—1.0 x 1078, + 1.0 x 107°] 
Hz/s (Abbott et al., 2017a). A separate Einstein@Home search using the GCT- 
F-statistic method drilled deeper in the 20-100 Hz band in a narrower spin- 
down range [—2.65 x 10~°,+2.64 x 10719] Hz/s (Abbott et al., 2017d). A 
follow-up publication using three of the first four pipelines (PowerFlux, Sky 
Hough and TD-F-statistic) covered the broader band 475-2000 Hz (Abbott 
et al., 2018b). 


Figure 31 shows the full-band O1 results from (Abbott et al., 2018b) for 
these three pipelines. The PowerFlux results shown are defined differently 
from those shown for the other searches. PowerFlux upper limits are derived 
as strict frequentist over the full sky, that is, a 95% CL limit provides at least 
95% coverage, regardless of sky position, making it quite conservative. At the 
same time, however, limits are shown for an optimistic polarization assump- 
tion (circular polarization corresponding to | cos(v)| = 1) and for a pessimistic 
assumption (linear polarization corresponding to cos(z) = 0 for the least favor- 
able choice of polarization angle ~). These limits are derived directly from the 
corresponding detection statistics (see section 3.5.2). The other limits shown 
are conventional frequentist population-based values, averaged over source ori- 
entation and sky position. Figure 32 shows the low-frequency band up to 100 
Hz, comparing the limits obtained in (Abbott et al., 2017a) with those from 
the GCT-F-statistic search on Einstein@Home, where the PowerFlux limits 
have been reevaluated via explicit simulation to produce population-averaged 
values for comparison. 
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All-sky results from three pipelines (Sky Hough, Frequency Hough and TD- 
F-statistic) were applied to the O2 data set (Abbott et al., 2019a; Palomba 
et al., 2019) over the 20-1922 Hz band and a spin-down range [—1.0x10~*, +2.0x 
10~°], where frequency coverage varied by pipeline. Resulting upper limits are 
shown in Figure 33. A dedicated Einstein@Home search of the O2 data (Stelt- 
ner et al., 2021) over the 20-585 Hz band achieved significantly lower upper 
limits in the overlapping frequency band (see Fig. 34). 

An intriguing set of Ol and O2 all-sky searches using the Falcon pipeline 
(derived from PowerF lux, but implemented with approximations and exploit- 
ing additional symmetries (Dergachev and Papa, 2019, 2020a,b, 2021a,b), fo- 
cused on deeper searches. The O1 search (Dergachev and Papa, 2019, 2020a) 
doubled the first-stage effective coherence time from that used in the O1 Pow- 
erFlux search (Abbott et al., 2017a, 2018b) while covering the same spin-down 
range over the 100-600 Hz band. The O2 searches, on the other hand, targeted 
low-ellipticity pulsars (Woan et al., 2018) by severely restricting the spin-down 
range (e.g., |faw| < 3 x 107!? Hz/s) in the 500-1500 Hz band). This vast re- 
duction in parameter space permits using loose coherence with an effective 
coherence time of 12 hours in its initial search stage, albeit with a necessarily 
reduced astrophysical range because of the restricted spin-down restriction. 

Another deep O2 search (Wette ct al., 2021) focused on the narrow 171-172 
Hz band while restricting spin-down magnitudes below ~3 x 107'8 Hz/s. This 
search used a semi-coherent F-statistic technique with Graphics Processing 
Unit acceleration in the F-statistic computation), where the frequency band 
chosen was meant to optimize probability density of detection in a narrow 
band based on detector sensitivity and the known pulsar population. 

The first all-sky search of O3 data for isolated CW sources (Abbott et al., 
2021b) used the PowerFlux pipeline to examine the O3a data for the same 
broad parameter space in frequency and spin-down as used in the O1 search. 
A comparison of upper limits obtained from several O2 searches with those 
obtained from O8a search are shown in Fig. 34. Figure 35 shows the corre- 
sponding parameter space coverages. 

The most sensitive all-sky results to date for broad coverage of both fre- 
quency and spin-down were obtained recently from the full O3 data from three 
pipelines (Sky Hough, Frequency Hough and TD-F-statistic) (Abbott et al., 
2022a) and are shown in Fig. 36, in comparison with the O3a PowerFlux re- 
sults (Abbott et al., 2021b) and with the results from the new Viterbi-based, 
less sensitive but blazing-fast, SOAP pipeline. Also shown are recent O3a Fal- 
con results (Dergachev and Papa, 2022) over a restricted spin-down range. 
Figure 37 shows a comparison of the parameter space coverages of these dif- 
ferent searches. 


4.5 All-sky searches for binary stars 


Several methods have been proposed and implemented for carrying out a CW 
all-sky binary search. The first method, which was used in a published search 
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Fig. 31 O1 all-sky upper limits (95% CL) on ho for isolated stars from three semi-coherent 
search pipelines over the band 20-2000 Hz (Abbott et al., 2018b). The limits shown for the 
PowerFlux method correspond to best-case (circular polarization) and worst-case (linear 
polarization) over the entire sky, while the limits shown for the time-domain F-statistic and 
SkyHough methods correspond to population averages over the sky and source orientations. 
The steps in sensitivty apparent in the limits correspond to reductions in FFT coherence 
time as frequencies increase. 


of initial LIGO S6 data (Aasi et al., 2014a) is known as TwoSpect (Goetz 
and Riles, 2011). The program carries out a semi-coherent search over an 
observation time long compared to the maximum orbital period considered, 
while using coherence times short with respect to the orbital period. Fourier 
transforms are carried out over each row (fixed frequency bin) in a ~year- 
long spectrogram, and the resulting frequency-frequency plot is searched for 
characteristic harmonic patterns. 

Another developed pipeline, known as Polynomial, (van der Putten et al., 
2010) searches coherently using matched filters over an observation time short 
compared to the minimum orbital period considered. A bank of frequency 
polynomials in time is used for creating the matched filters, where for a 
small segment of an orbit, the frequency should vary as a low-order polyno- 
mial. Other proposed methods, which offer potentially substantial computa- 
tional savings at a cost in sensitivity, include autocorrelations in the time- 
frequency plane (Viceré and Yvert, 2016) and stochastic-background tech- 


116 Keith Riles 


4 Powerflux O1 search 

= Time-domain F-stat O1 search 
» Sky Hough O1 search 
Frequency Hough O1 search 
Results from this search 


10-74 


90% 
hg 


10-25 


20 30 40 50 60 70 80 90 100 
search frequency (Hz) 


Fig. 32 O1 all-sky upper limits (95% CL) on ho for isolated stars in the low-frequency 
band (20-100 Hz) for five semi-coherent pipelines (Abbott et al., 2017d), including an Ein- 
stein@Home GCT-F-statistic search. The PowerFlux limits here are population-averaged, 
unlike those shown in Figure 31. 


niques (Ballmer, 2006b), with computational costs gains achieved by using 
skymaps with sidereal-day folding (Thrane et al., 2015; Goncharov and Thrane, 
2018; Abbott et al., 2021a). 


More recently, the implementation of graphics processor units software in 
the framework of the Sky Hough all-sky program has led to a breakthrough in 
all-sky binary search sensitivity (Covas and Sintes, 2019). Upper limits were 
initially obtained over 100-300 Hz and over a broad range of binary orbital 
parameters from the LIGO O02 data(Covas and Sintes, 2020). Although this 
approach does not yet cover the full orbital parameter space possible with the 
TwoSpect program, the intrinsic sensitivity is dramatically better, with exten- 
sion of the method to shorter orbital periods a natural future improvement. 
A follow-up analysis in the O3a data (Abbott et al., 2021c; Tenorio, 2021) 
expanded the search band slightly (50-300 Hz), and a parallel development 
using a similar Hough transform framework but with a F-statistic (Covas and 
Prix, 2022) tailored to multi-hour segments, has been applied to the O3a data 
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Fig. 33 O2 all-sky upper limits (95% CL) on ho for isolated stars from three semi-coherent 
search pipelines over the band 20-1922 Hz. As in Figure 31, step changes in sensitivity 
correspond to reductions in FFT coherence time with increasing frequency. 


in the 300-500 Hz band (Covas et al., 2022). All of these results are shown in 
Figure 38. See references for details on orbital parameter space regions covered 
by the different analyses, which vary considerably. 

In addition, searches for isolated stars retain some sensitivity to long-period 
binaries, as detailed in a recent study (Singh et al., 2019). 


4.6 Searches for CW transients and other CW-like signals 


The first dedicated search for CW transients following a known pulsar’s glitch 
addressed glitches detected by radio astronomers during the Advanced LIGO 
O2 data run. The search used the transient F-statistic method (Prix et al., 
2011) and focused on periods following glitches in the Crab and Vela pul- 
sars (Keitel et al., 2019). A recent 03 analysis (Abbott et al., 2021h; Modaf- 
feri et al., 2021) searched for CW transients following nine glitches across 
six pulsars (one glitch each from five stars: PSR J0534+2200, J0908-4913, 
J1105-6107, J1813-1749 and J1826-—1334; and four glitches from the intrigu- 
ing source, PSR, J0537-6910 (see sections 2.1.2, 2.1.4 and 4.2). No significant 
candidates were observed, although two marginal outliers were seen after one 
PSR J0537-6910 glitch, albeit with implied strengths well above those con- 
sistent with the inferred glitch energies. In fact, the upper limits obtained for 
post-glitch energy emission from all glitches examined lay above the maximum 
expected in a simple two-fluid model (Prix et al., 2011), with strain limits for 
PSR J1105-6107 approaching most closely to that benchmark (within a factor 
of ~1.6). 
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Fig. 34 O3a all-sky upper limits (95% CL) on ho for isolated stars from the 03a PowerF lux 
search in comparison with earlier O2 searches. The corresponding parameter space areas in 
faw-—faw are shown in Fig. 35. 


The first dedicated search for long-lived, CW-like signals from a post- 
merger remnant looked for a signal from the post-GW170817 remnant, but as 
expected, given the ~40 Mpc distance to the merger, no signal was detected in 
the immediate aftermath (Abbott et al., 20170) (~500 s) or in a multi-hour to 
multi-day period afterward (Abbott et al., 2019d). Should another opportu- 
nity arise (from a nearby binary neutron star merger or a galactic supernova), 
search methods are available for use (Thrane et al., 2011; Miller et al., 2018; 
Sun and Melatos, 2019; Oliver et al., 2019; Banagiri et al., 2019; Mytidis et al., 
2019; Miller et al., 2019a). 

More exotic recent analyses seeking CW or CW-like signals include: 


— Searches for CW signals from Bose-Einstein clouds (D’ Antonio et al., 2018; 
Abbott et al., 2022b) (see section 2.2). 

— A search for non-black-hole weakly interacting compact dark objects with 
mass below 10~" Mo orbiting within the Sun about its center (Horowitz 
et al., 2020) 

— Searches for ultralight dark photon or scalar boson dark matter creating 
an extremely narrowband (Af/f ~ 10~°) spectral excess with stochastic 
phase (Pierce et al., 2018; Guo et al., 2019; Miller et al., 2021b; Abbott 
et al., 2022c; Grote and Stadnik, 2019; Vermeulen et al., 2021). 
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Fig. 35 Comparison of parameter space areas in O2 all-sky searches vs the O3a PowerFlux 
search. The shaded rectangle with vertical bars shows the 20-2000 Hz and —10~8-10-° Hz/s 
range for the O3a search (Abbott et al., 2021b). The slightly larger rectangle with horizontal 
bars shows the region searched in the O2 data with the Frequency Hough method (Abbott 
et al., 2019a; Palomba et al., 2019). The smaller rectangle with crossed diagonal bars shows 
the region searched by the distributed-computing project Einstein@Home (Stceltner et al., 
2021). The solid line at zero spin-down depicts the specialized O2 search for low-ellipticity 
millisecond pulsars using the Falcon method (Dergachey and Papa, 2020b, 2021a,b) (the 
thickness of the line overstates the coverage in spin-down range). The dotted curves indicate 
contours of constant equatorial ellipticity « = (10-8, 10-7, 10—®, 10-5 and 10~*) for a star 
with stellar spin-down dominated by gravitational wave emission. 


— A search for binary systems of planetary-scale / asteroid-scale primordial 
black holes near to the Earth (Miller et al., 202la, 2022; Abbott et al., 
2022a). 

— There has also been a proposal to apply CW search techniques to suspected 
Thorne-Zytkow objects (Thorne and Zytkow, 1975) (TZOs) for which a 
neutron star orbiting inside of a giant star (slow inspiral decay) could pro- 
duce signals in the band of ground-based gravitational wave detectors (De- 
Marchi et al., 2021). 


5 Outlook 
5.1 Prospects for discovery 


Over the next several years, the Advanced LIGO (Aasi et al., 2015a) and 
Virgo (Acernese et al., 2014) detectors are expected to approach and eventu- 
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Fig. 36 O3 all-sky upper limits (95% CL) on ho for isolated stars from four pipelines (Ab- 
bott et al., 2022a), in comparison with the O3a PowerFlux (Abbott et al., 2021b) and Falcon 
results (Dergachev and Papa, 2022). The corresponding parameter space areas in faw-faw 
are shown in Fig. 37. 


ally surpass their original design sensitivities in strain, increasing the range 
within the galaxy which CW searches can access, thereby increasing detection 
likelihood. As the sensitive ranges of different search methods approach the 
dense galactic core, detection chances may rise more rapidly. In parallel to de- 
tector improvement, algorithms continue to improve, as researchers find more 
effective tradeoffs between computational cost and detection efficiency, while 
Moore’s Law, including Graphics Processing Unit (GPU) exploitation (Keitel 
and Ashton, 2018; Covas and Sintes, 2019; Wette et al., 2021; Dunn et al., 
2021; Rosa et al., 2021), ensures increased computing resources for searches. 
All of these trends are encouraging for successful CW detection. 

At the same time, theoretical uncertainties in what sensitivity is needed 
for the first CW detection are very large. While the spin-down limits based 
on gravitar assumptions and on either energy conservation or known age have 
been beaten for a handful of sources and will be beaten for more sources in 
the coming years, the gravitar model is surely optimistic — most stellar spin- 
downs are likely dominated by electromagnetic interactions. Whether the first 
detection is imminent or still many years distant remains unclear. A recent 
phenomenological population synthesis study (Cieslar et al., 2021), based on 
an exponentially decaying ellipticity that starts at its allowed maximum ~107° 
with a supernova rate of once per century concluded that the expected number 
of detectable, young and isolated neutron stars for Advanced LIGO sensitivity 
is less than one and is ~10 for Einstein Telescope. 

Electromagnetic astronomers could prove pivotal in hastening detection 
by identifying new nearby or young neutron stars, or discovering pulsations 
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Fig. 37 Comparison of faw-—faw parameter space coverage for the four search pipelines 
used in the full-O3 all-sky searches (Abbott et al., 2022a) and for the restricted-spindown 
Oa Falcon search. (Abbott et al., 2022a). 


from known stars, perhaps most usefully from the accreting Sco X-1 sys- 
tem (Galaudage et al., 2021). Given the computational challenges of most 
CW searches, narrowing the parameter space of a search exploiting electro- 
magnetic observations could make the difference between a gravitational wave 
miss and a discovery. 


5.2 Confirming and exploiting a discovery 


There are several aspects of confirming a nominal CW discovery, including 
establishing the statistical significance of the outlier, verifying consistency of 
the signal with the CW model, excluding an environmental or instrumental 
cause, and (optionally but ideally) confirming consistency with prior or follow- 
up electromagnetic observations. Once that discovery is established, exploiting 
it to understand neutron star astrophysics (or fundamental particle physics 
should a superradiant boson cloud be observed) will be a rich endeavor. 

The degree of statistical confidence with which a putative CW signal de- 
tection can be confirmed depends on the type of search that leads to the 
candidate. The statistical significance of an outlier depends on a trials factor 
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Fig. 38 All-sky upper limits (95% CL) on ho for stars in binary systems. Upper limits 
are shown from the inital LIGO S6 TwoSpect search (Aasi et al., 2014a), from the GPU- 
enhanced O2 Binary Sky Hough search (Covas and Sintes, 2020) (100-300 Hz), from the 
O3a Binary Sky Hough search (Abbott et al., 2021c) (50-300 Hz) and from the O3a Binary 
F-statistic search (Covas ct al., 2022) (300-500 Hz). See references for details on orbital 
parameter space regions covered by the different analyses. The O2 and O3a Binary Sky 
Hough values shown are 95% sensitivities with bands to indicate uncertainties. 


that may range from ~1 for targeted searches for known pulsars using elec- 
tromagnetically derived ephemerides to ~10!° for all-sky searches. Hence the 
SNR threshold for, say, a “5-sigma” discovery varies too. In practice, though, 
for a candidate emerging from a hierarchical search with multiple stages of 
“zooming in”, the SNR for a surviving outlier may be so much higher than the 
sensitivity-defining SNR threshold used in the first stage that the initial trials 
factor is irrelevant. Establishing the statistical confidence of a targeted-search 
candidate, on the other hand, may simply require steady accumulation of ad- 
ditional data while fully exploiting all data in hand from all detectors with 
appreciable sensitivity in that band. Empirically assessing the significance of 
loudest outliers is discussed in detail in (Tenorio et al., 2022). 


If known detector artifacts are degrading sensitivity in the frequency band 
of the candidate, it may be feasible to focus detector commissioning to miti- 
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gate the artifact prior to the next observing run. Another possible approach, 
although potentially detrimental to other GW observations, is “narrowband- 
ing.” In narrowbanding the detector sensitivity is improved in a narrow band 
at the expense of broadband sensitivity by adjusting the position of the “signal 
recycling” mirror at the output port of the interferometer (Aasi ct al., 2015a); 
with the advent of quantum squeezing in advanced detectors (Tse et al., 2019), 
however, the potential gains from narrowbanding are less pronounced. 

Even with a promising outlier, a discovery claim would need more than sta- 
tistical inconsistency with detector noise. One would seek consistency with the 
signal model, particularly for candidates originating in hierarchical searches 
where early stages look primarily for excess power that is only roughly consis- 
tent with a particular template and where the spacing between templates is 
relatively coarse. The expected Doppler modulations due to the Earth’s mo- 
tion should be present (Zhu et al., 2017; Intini et al., 2020b). One wants to see 
a signal for which a fully coherent search over all data yields an SNR consis- 
tent with expectation from the putative source. Ideally, a residual spectrum 
from subtracting the reconstructed signal would be consistent with random 
background noise. 

In the case of a targeted or directed search for which the source location 
is a priori known, one would want to verify that the highest-SNR template 
observed in that region of the sky and near the template’s frequency param- 
eters is indeed consistent with the correct sky location. Although one could 
impose a similar constraint on frequency and frequency derivatives for a tar- 
geted search, narrowband searches do allow those parameters to differ slightly 
from the nominal ones, to accommodate differential stellar rotation. Hence 
seeing the SNR peak at precisely the right location in parameter space for a 
known pulsar would lend credence to the signal, but seeing the SNR peak at a 
nearby point in parameter space can still mean discovery, albeit with a trials 
factor appropriate to a narrowband search. 

The possibility of a rotational glitch during an observational period presents 
additional complications. One can no longer safely apply a fully coherent search 
over the time span and expect a monotonically increasing SNR. In the case of 
a targeted search with ephemerides in hand indicating a glitch, breaking the 
observation time into two (or more) segments is straightforward (Abbott et al., 
2010), but in the case of a source without independent timing information, one 
may have a true signal detection but lack the confidence to declare discovery 
without additional data taking because of apparent phase inconsistency in the 
available data. 

In any gravitational wave analysis using interferometers that push the fron- 
tier of technology (and which are routinely operated at maximum achievable 
sensitivity), one must consider whether or not instrumental or environmen- 
tal contamination leads to a false signal. As discussed in section 3.7, narrow 
lines can contribute to accumulated power in a templated search. As part of 
confirming a discovery, one would need to quantify that contamination for a 
putative signal lying near a known instrumental spectral line. More challenging 
and more realistic for a signal candidate surviving multiple hierarchical search 
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stages are spectral lines that are not immediately apparent in the strain chan- 
nel spectrum, especially lines that are non-stationary with respect to time or 
frequency (“wandering”). To address that possibility, one would look compre- 
hensively at auxiliary data channels, such as readouts from magnetometers, 
accelerometers, seismometers, microphones and from any servo control chan- 
nels that could impose tiny actuations on the gravitational wave strain chan- 
nel. Those investigations would include examination of averaged spectra for 
peaks coincident with the strain signal frequency and more probing searches 
for cross-coherence between the auxiliary channels and the strain channel that 
is inconsistent with statistical fluctuation. 

Finally, in confirming a continuous gravitational wave signal one would, 
ideally, want confirmation via electromagnetic observations. For targeted or 
narrowband searches of known pulsars, the observations already exist, and 
the primary task is to establish statistical confidence of their consistency with 
gravitational data. For other known sources, however, gravitational wave mea- 
surements may provide the necessary clues to allow detection of previously un- 
detected pulsations. For example, detection of a CW signal from Scorpius X-1 
could permit discrimination of X-ray pulsations from a stochastic background 
dominated by accretion emission. For a previously unknown source found in 
an all-sky search, determining the source location from coherent integration 
over months of data (with sub-arcsecond resolution possible from the aperture 
formed by the Earth’s orbit) may suffice for radio, X-ray, gamma-ray (and 
perhaps even optical) astronomers to find the counterpart. If electromagnetic 
pulsations were detected and agreed with expectation, the confirmation of the 
gravitational wave signal would be ironclad. 

An interesting challenge to confirmation would be continuous gravitational 
radiation due to boson cloud superradiance for an isolated black hole (see sec- 
tion 2.2). If there were no accretion disk or companion to induce an electro- 
magnetic signal, one would have to rely heavily upon the the evolution of the 
gravitational wave signal itself to infer the nature of the source. In particular, 
the source frequency governed by the boson’s apparent mass in the potential 
of the black hole could spin up instead of down as the black hole loses mass 
energy to gravitational radiation, thereby reducing the magnitude of the neg- 
ative binding energy correction to the unbound boson mass (Arvanitaki et al., 
2015). 

Once a continuous gravitational wave detection has been confirmed elec- 
tromagnetically, one will want to exploit the correlations to understand the 
source. Below are a sampling of potential measurements possible, along with 
questions they help to address: 


— Relation between rotational and gravitational wave frequencies, determin- 
ing the fundamental mechanism of emission (see section 2 and see (Jones, 
2021) for a detailed discussion). 

— Correlation of the gravitational and electromagnetic phase constants in 
the event of consistent frequency (phase) evolution. If an equatorial mass 
“bulge” explains the GW signal, how well does the implied quadrupolar axis 
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align with a pulsar’s inferred magnetic dipole projection? In an accreting 
system, for example, does added mass accumulate near the magnetic poles? 

— Differential frequency (phase) evolution. Is there differential rotation be- 
tween the stellar crust and its interior? If electromagnetic frequency glitches 
are observed, what is seen gravitationally before, during and after the 
glitch? 

— If there is evidence for r-modes from, say, an approximate 4/3 ratio of 
GW signal frequency to stellar rotation frequency, how does the GW fre- 
quency evolve with time and how does the ratio evolve? Is there evidence 
of amplitude growth from instability? Decay from viscosity? 

— Inferred quadrupole moment. Although ellipticity is a convenient dimen- 
sionless parameter, it is approximately the product of the ellipticity times 
the stellar moment of inertia about its spin axis that determines the signal 
strength for a mass-quadrupole radiator. Given the uncertainties in neu- 
tron star equation of state, there are large uncertainties in the moment of 
inertia and hence ambiguity in extracting ellipticity. Ambiguity at the level 
of near-degeneracy would arise in the absence of an independent determi- 
nation of source distance from electromagnetic observations (Sicniawska 
and Jones, 2021). A CW detection in a binary system would offer an op- 
portunity for determination of the stellar mass. Other stellar properties po- 
tentially accessible include the stellar radius (inferred from luminosity and 
temperature, if measurable). A precessing star with detectable electromag- 
netic pulsations offers additional opportunity for understanding internal 
structure (Gao et al., 2020). 

— Boson properties from superradiance. In the event of detecting superradi- 
ance from a boson cloud around a black hole, determining the boson’s mass 
will be immediate from the signal frequency (at least for the annihilation 
channel expected to dominate) with the boson intrinsic spin determination 
more model dependent, based on signal strength and frequency evolution 
with some knowledge of the black hole source needed. 


In addition to exploiting CW detection to understand the source, one can 
also carry out precise tests of General Relativity by measuring the polarization 
of the propagating gravitational wave. In Einstein’s theory there should be 
two independent transverse, quadrupolar polarizations for which the relative 
strengths depend on the source orientation relative to the line of sight. In non- 
standard theories of gravity other polarization modes, including scalar, vector 
and longitudinal polarizations, may be present (Isi ct al., 2015). Testing for 
these additional polarizations with transient gravitational wave detections to 
date has been challenging because nearly all of the signal-to-noise ratio has 
come from the two nearly aligned LIGO detectors such that they mainly detect 
the same polarization projection. In contrast, even for a single detector, a CW 
signal permits disentangling multiple polarization contributions as the sidereal 
rotation of the Earth changes the detector’s (polarimeter’s) orientation with 
respect to the source direction deterministically (Isi et al., 2015; Kuwahara 
and Asada, 2022). In fact, even in the absence of a CW signal, one can set 
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upper limits on the non-standard polarizations (Abbott et al., 2018a), just as 
is possible for standard polarizations. 

Nature has blessed the gravitational wave community with a bounty of 
compact star mergers, including the remarkable first detected BBH merger, 
GW150914, and the even more remarkable and informative multi-messenger 
detection of the GW170817 BNS event. Should such kindness continue, one 
may hope soon for a multi-messenger detection of a CW source that not only 
could be observed into the foreseeable future, but could mark the first of a 
large collection to come, as GW150914 proved to be. 
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